1 Setup

2 Import

data <- read.csv2("data/d.csv", header = TRUE)

3 Analysis

Just a short information on the amount of pause values to be found in the ungrouped condition to find out, if it is reasonable to model pause data for this condition.

# how many not zero pauses do we have in the ungrouped condition  ---------------------------------
grou <- subset(data, data$condition=="grouped")
colSums(grou != 0)
##          X.1            X         Subj        group          exp          rep 
##          658          658          658          658          658          658 
##        trial    condition    f0_range1    f0_range2      s4reln1      s8reln2 
##           NA          658           NA           NA          658          658 
##       pause1       pause3       meanf0         list           L1          L1t 
##           80          522          658          658           NA           NA 
##           H1          H1t         f0n1      slopen1           L2          L2t 
##           NA           NA           NA           NA           NA           NA 
##           H2          H2t         f0n2      slopen2        s4dur        s8dur 
##           NA           NA           NA           NA          658          658 
##        p1dur        p2dur        p3dur        p4dur       fp0dur       fp1dur 
##           80           53          522          141            1            0 
##       fp2dur       fp3dur       fp4dur        c1dur        c2dur  utt_dur.ms. 
##            7            0            3          658          658          658 
##     filename rating_match         item       pause2       pause4     accuracy 
##          658           NA          658           53          141          658 
##         rate 
##          428
rm(grou)
un_grou <- subset(data, data$condition=="ungrouped")
colSums(un_grou != 0)
##          X.1            X         Subj        group          exp          rep 
##          643          643          643          643          643          643 
##        trial    condition    f0_range1    f0_range2      s4reln1      s8reln2 
##           NA          643           NA           NA          643          643 
##       pause1       pause3       meanf0         list           L1          L1t 
##          118          186          643          643           NA           NA 
##           H1          H1t         f0n1      slopen1           L2          L2t 
##           NA           NA           NA           NA           NA           NA 
##           H2          H2t         f0n2      slopen2        s4dur        s8dur 
##           NA           NA           NA           NA          643          643 
##        p1dur        p2dur        p3dur        p4dur       fp0dur       fp1dur 
##          118           38          186           65            1            0 
##       fp2dur       fp3dur       fp4dur        c1dur        c2dur  utt_dur.ms. 
##            6            1            4          643          643          643 
##     filename rating_match         item       pause2       pause4     accuracy 
##          643           NA          643           38           65          643 
##         rate 
##          509
rm(un_grou)
  • non-zero pause values in ungrouped condition:
    • data points overall: 643
      • pause1: 118
      • pause2: 38
      • pause3: 186 (29% of data points overall in ungrouped condition)
      • pause4: 65
  • non-zero pause values in grouped condition:
    • data points overall: 658
      • pause1: 80
      • pause2: 53
      • pause3: 522 (79% of data points overall in grouped condition)
      • pause4: 141

I would suggest, it makes sense to include pause3 in both conditions.

3.1 Contrast coding

  • assign scaled sum contrast to all factors - all have two levels (0.5, -0.5)
  • intercept estimates average effects
  • slopes estimate difference between factor levels
# Contrast coding  ---------------------------------
## we have 4 predictors / factors with two levels (condition, group, exp, accuracy)
## therefore it is best to use sum contrast with 0.5 / -0.5
## within each predictor this compares averages of each level against each other - mean of these contrasts is 0
## interactions between predictors estimate average effects (estimates are at mean which is 0 = averaged)
# Prepare data  ---------------------------------
## Subset for patient data - only RHDP / LHDP =================================
## model stimuli (stimuli that had to be repeated) will be "removed"
pDat <- subset(data, data$group!="model_stimuli")
### Factorize predictors variables #############################
########################################### condition
class(pDat$condition) #character
## [1] "character"
pDat$condition <- as.factor(pDat$condition)
## check level order 
levels(pDat$condition) # "grouped" "ungrouped"
## [1] "grouped"   "ungrouped"
########################################### group
class(pDat$group) #character 
## [1] "character"
pDat$group <- as.factor(pDat$group)
## check level order 
levels(pDat$group) # "LHDP", "RHDP" 
## [1] "LHDP" "RHDP"
########################################### exp
class(pDat$exp) # character
## [1] "character"
pDat$exp <- as.factor(pDat$exp)
## check level order
levels(pDat$exp) # "reading_aloud" "repetition" 
## [1] "reading_aloud" "repetition"
########################################### accuracy
class(pDat$accuracy) # "character"
## [1] "character"
## factorize "accuracy"
pDat$accuracy <- as.factor(as.character(pDat$accuracy))
## check level order
levels(pDat$accuracy) # "correct"   "incorrect"
## [1] "correct"   "incorrect"
### contrast code predictor variables #############################
########################################### condition
contrasts(pDat$condition) <- contr.sum(2)/2
contrasts(pDat$condition)
##           [,1]
## grouped    0.5
## ungrouped -0.5
#             [,1]
# grouped     0.5
# ungrouped -0.5
########################################### group
## define sum contrast for group
contrasts(pDat$group) <- contr.sum(2)/2
contrasts(pDat$group)
##      [,1]
## LHDP  0.5
## RHDP -0.5
#       [,1]
# LHDP  0.5
# RHDP -0.5
########################################### exp
## assign sum contrast to experiment
contrasts(pDat$exp) <- contr.sum(2)/2
contrasts(pDat$exp)
##               [,1]
## reading_aloud  0.5
## repetition    -0.5
#             [,1]
# reading_aloud      0.5
# repetition  -0.5
########################################### accuracy
## define sum contrast for accuracy
contrasts(pDat$accuracy) <- contr.sum(2)/2
contrasts(pDat$accuracy)
##           [,1]
## correct    0.5
## incorrect -0.5
#           [,1]
# correct    0.5
# incorrect -0.5

3.2 Find correct / interpretable / parsimonious model to analyze data

3.2.1 Maximally complex model based on our hyoptheses

  • the random effects structure is reduced iteratively in the next step to avoid overparameterization as best as possible
# Max models based on hypotheses ---------------------------------
## FL name2 =================================
summary(FL2_h <- lmer(s8reln2 ~ condition*accuracy+condition*group*exp+
              group*exp*accuracy+
              (1+condition*accuracy+exp*accuracy+condition*exp|Subj)+
              (1+group*accuracy|item),
            REML=TRUE,pDat))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: s8reln2 ~ condition * accuracy + condition * group * exp + group *  
##     exp * accuracy + (1 + condition * accuracy + exp * accuracy +  
##     condition * exp | Subj) + (1 + group * accuracy | item)
##    Data: pDat
## 
## REML criterion at convergence: 8222.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2221 -0.6159 -0.0027  0.5753  4.0519 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr                         
##  Subj     (Intercept)          18.894   4.347                                 
##           condition1            7.272   2.697    -0.72                        
##           accuracy1             2.962   1.721    -0.47  0.29                  
##           exp1                 10.784   3.284     0.05 -0.07  0.47            
##           condition1:accuracy1 51.198   7.155    -0.14  0.18 -0.80 -0.43      
##           accuracy1:exp1       13.125   3.623    -0.02  0.19 -0.03 -0.10  0.05
##           condition1:exp1      12.211   3.494    -0.13 -0.19 -0.30 -0.24  0.39
##  item     (Intercept)           7.809   2.794                                 
##           group1                5.637   2.374    -0.99                        
##           accuracy1             3.373   1.837     0.83 -0.88                  
##           group1:accuracy1     10.218   3.197     0.92 -0.96  0.98            
##  Residual                      27.753   5.268                                 
##       
##       
##       
##       
##       
##       
##       
##  -0.84
##       
##       
##       
##       
##       
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                        Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)             43.4670     0.9944 45.3883  43.711  < 2e-16 ***
## condition1               4.7019     1.2593 31.3445   3.734 0.000752 ***
## accuracy1               -0.1099     0.6756 28.6661  -0.163 0.871936    
## group1                   3.9785     1.6936 33.0783   2.349 0.024943 *  
## exp1                     1.7234     1.3869 31.5163   1.243 0.223160    
## condition1:accuracy1     9.5060     1.6397 26.0327   5.797 4.14e-06 ***
## condition1:group1       -0.4442     1.5411 38.2497  -0.288 0.774703    
## condition1:exp1          0.2030     2.2923 29.2689   0.089 0.930022    
## group1:exp1             -4.2497     1.7859 35.0612  -2.380 0.022904 *  
## accuracy1:group1        -1.7295     1.2535 35.4761  -1.380 0.176306    
## accuracy1:exp1          -0.7491     1.3584 26.8565  -0.551 0.585889    
## condition1:group1:exp1  -2.2148     2.1714 31.6666  -1.020 0.315482    
## accuracy1:group1:exp1    4.3390     2.6116 29.5572   1.661 0.107204    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# Model failed to converge with max|grad| = 0.0250477 (tol = 0.002, component 1)
# reduction needed
## F0 name2 =================================
summary(F02_h <- lmer(f0_range2 ~ condition*accuracy+condition*group*exp+
              group*exp*accuracy+
              (1+condition*accuracy+exp*accuracy+condition*exp|Subj)+
              (1+group*accuracy|item),
            REML=TRUE,pDat))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: f0_range2 ~ condition * accuracy + condition * group * exp +  
##     group * exp * accuracy + (1 + condition * accuracy + exp *  
##     accuracy + condition * exp | Subj) + (1 + group * accuracy |      item)
##    Data: pDat
## 
## REML criterion at convergence: 5846.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4215 -0.4556 -0.0352  0.3921  6.9327 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr                         
##  Subj     (Intercept)          3.32418  1.8232                                
##           condition1           1.75620  1.3252    0.38                        
##           accuracy1            0.22160  0.4707    0.44  1.00                  
##           exp1                 1.84342  1.3577   -0.14 -0.02 -0.03            
##           condition1:accuracy1 6.70690  2.5898    0.54 -0.06 -0.03 -0.38      
##           accuracy1:exp1       1.25981  1.1224    0.72  0.64  0.67 -0.14  0.17
##           condition1:exp1      2.33252  1.5273   -0.03  0.56  0.54  0.77 -0.57
##  item     (Intercept)          0.04047  0.2012                                
##           group1               0.62260  0.7890   -0.82                        
##           accuracy1            0.16225  0.4028    0.32 -0.80                  
##           group1:accuracy1     0.45821  0.6769    0.95 -0.59  0.00            
##  Residual                      5.42452  2.3291                                
##       
##       
##       
##       
##       
##       
##       
##   0.24
##       
##       
##       
##       
##       
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                        Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)              5.2720     0.3446 29.9811  15.298 1.05e-15 ***
## condition1               2.0080     0.3113 27.0110   6.451 6.50e-07 ***
## accuracy1                0.3662     0.2373 52.0728   1.543 0.128914    
## group1                   0.9976     0.6488 36.0793   1.538 0.132854    
## exp1                     2.2673     0.3369 27.3297   6.730 2.98e-07 ***
## condition1:accuracy1     3.0841     0.6293 27.2214   4.900 3.90e-05 ***
## condition1:group1        0.1046     0.6575 31.1750   0.159 0.874670    
## condition1:exp1          0.8512     0.4699 34.8730   1.811 0.078710 .  
## group1:exp1              2.8110     0.7055 34.3688   3.984 0.000334 ***
## accuracy1:group1        -0.2812     0.4659 44.2023  -0.604 0.549238    
## accuracy1:exp1           0.3828     0.4848 25.3706   0.790 0.437107    
## condition1:group1:exp1   1.8709     0.9410 36.0065   1.988 0.054433 .  
## accuracy1:group1:exp1    1.4779     0.9517 21.8607   1.553 0.134785    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
# reduction needed
## Pause3 =================================
summary(P_h <- lmer(pause3 ~ condition*accuracy+condition*group*exp+
            group*exp*accuracy+
            (1+condition*accuracy+exp*accuracy+condition*exp|Subj)+
            (1+group*accuracy|item),REML=TRUE,pDat))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pause3 ~ condition * accuracy + condition * group * exp + group *  
##     exp * accuracy + (1 + condition * accuracy + exp * accuracy +  
##     condition * exp | Subj) + (1 + group * accuracy | item)
##    Data: pDat
## 
## REML criterion at convergence: 7906.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4980 -0.3586 -0.0742  0.3440  6.8107 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr                         
##  Subj     (Intercept)           7.89696 2.8102                                
##           condition1            9.62927 3.1031    0.24                        
##           accuracy1             5.12554 2.2640   -0.20  0.43                  
##           exp1                 14.10797 3.7561    0.27  0.38 -0.19            
##           condition1:accuracy1 54.71130 7.3967    0.78  0.07  0.08 -0.31      
##           accuracy1:exp1       11.64825 3.4130    0.33  0.33  0.70 -0.47  0.74
##           condition1:exp1      31.55293 5.6172    0.32  0.74  0.00  0.72  0.03
##  item     (Intercept)           0.00000 0.0000                                
##           group1                0.34719 0.5892     NaN                        
##           accuracy1             0.06886 0.2624     NaN  1.00                  
##           group1:accuracy1      6.89504 2.6258     NaN -1.00 -1.00            
##  Residual                      22.89651 4.7850                                
##       
##       
##       
##       
##       
##       
##       
##  -0.02
##       
##       
##       
##       
##       
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                        Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)              8.9994     0.5498 28.5552  16.370 4.78e-16 ***
## condition1               8.4077     0.6865 29.6301  12.246 4.03e-13 ***
## accuracy1                2.0065     0.6163 27.1148   3.256  0.00303 ** 
## group1                   2.0490     0.9164 31.5150   2.236  0.03257 *  
## exp1                     4.6392     0.8287 24.9723   5.598 8.03e-06 ***
## condition1:accuracy1    17.5473     1.5823 28.6336  11.090 7.05e-12 ***
## condition1:group1        1.1401     1.3838 31.1602   0.824  0.41628    
## condition1:exp1          7.1036     1.2255 25.4400   5.797 4.52e-06 ***
## group1:exp1             -0.7602     1.5975 28.2742  -0.476  0.63781    
## accuracy1:group1        -0.2328     1.3381 28.6351  -0.174  0.86309    
## accuracy1:exp1           0.3461     1.1330 31.3823   0.306  0.76199    
## condition1:group1:exp1  -0.8169     2.4463 26.0634  -0.334  0.74109    
## accuracy1:group1:exp1    2.0963     2.3228 32.2679   0.902  0.37348    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
# reduction needed

4 rePCA

4.1 Indicator variables for all parameters

Create model.matrix dependent on model specification to merge these cols for all indicator variables with data.frame to model data

# generate model matrix  ---------------------------------
pDat_mm <- pDat
mm <- model.matrix(~condition*accuracy+condition*group*exp+group*exp*accuracy, pDat)
colnames(mm)
##  [1] "(Intercept)"            "condition1"             "accuracy1"             
##  [4] "group1"                 "exp1"                   "condition1:accuracy1"  
##  [7] "condition1:group1"      "condition1:exp1"        "group1:exp1"           
## [10] "accuracy1:group1"       "accuracy1:exp1"         "condition1:group1:exp1"
## [13] "accuracy1:group1:exp1"
## create new col in pDat from mm cols =================================
# I renamed the col to reflect the level comparisons
## condition1 - grpd.ungr
## accuracy1 - cor.inc
## group1 - lhd.rhd
## exp1 - rd.rpt
pDat_mm$grpd.ungr                 <-mm[,2]
pDat_mm$cor.inc                   <-mm[,3]
pDat_mm$lhd.rhd                   <-mm[,4]
pDat_mm$rd.rpt                    <-mm[,5]
pDat_mm$grpd.ungr_cor.inc         <-mm[,6]
pDat_mm$grpd.ungr_lhd.rhd         <-mm[,7]
pDat_mm$grpd.ungr_rd.rpt          <-mm[,8]
pDat_mm$lhd.rhd_rd.rpt            <-mm[,9]
pDat_mm$cor.inc_lhd.rhd           <-mm[,10]
pDat_mm$cor.inc_rd.rpt            <-mm[,11]
pDat_mm$grpd.ungr_lhd.rhd_rd.rpt  <-mm[,12]
pDat_mm$cor.inc_lhd.rhd_rd.rpt    <-mm[,13]
# Clean workspace ---------------------------------
rm(mm)

4.2 Model reduction - final lengthening (FL2)

# all "REML = FALSE" to use ANOVA model comparison later on
# INITIAL model with indicator variables ---------------------------------
summary(m0_FL2 <- lmer(s8reln2 ~ 
                         1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                         grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                         cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                         lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                         grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                         (1+cor.inc+grpd.ungr+rd.rpt+
                            grpd.ungr_cor.inc+
                            cor.inc_rd.rpt+
                            grpd.ungr_rd.rpt|Subj)+
                         (1+lhd.rhd+cor.inc+
                            cor.inc_lhd.rhd|item),
                       REML = FALSE, pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## s8reln2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + cor.inc_rd.rpt + grpd.ungr_rd.rpt | Subj) +  
##     (1 + lhd.rhd + cor.inc + cor.inc_lhd.rhd | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8358.2   8626.6  -4127.1   8254.2     1237 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2136 -0.6172 -0.0114  0.5835  4.0524 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)       17.629   4.199                                 
##           cor.inc            2.609   1.615    -0.45                        
##           grpd.ungr          6.689   2.586    -0.73  0.26                  
##           rd.rpt             9.973   3.158     0.05  0.49 -0.06            
##           grpd.ungr_cor.inc 49.570   7.041    -0.14 -0.81  0.19 -0.44      
##           cor.inc_rd.rpt    10.508   3.242     0.01 -0.03  0.21 -0.06  0.03
##           grpd.ungr_rd.rpt  10.963   3.311    -0.14 -0.30 -0.20 -0.27  0.40
##  item     (Intercept)        6.752   2.598                                 
##           lhd.rhd            4.961   2.227    -1.00                        
##           cor.inc            2.836   1.684     0.88 -0.92                  
##           cor.inc_lhd.rhd    8.860   2.977     0.95 -0.97  0.99            
##  Residual                   27.763   5.269                                 
##       
##       
##       
##       
##       
##       
##       
##  -0.86
##       
##       
##       
##       
##       
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)               43.4630     0.9498 48.5381  45.761  < 2e-16 ***
## grpd.ungr                  4.6890     1.1872 35.8285   3.950 0.000351 ***
## lhd.rhd                    3.9835     1.6340 35.0837   2.438 0.019978 *  
## cor.inc                   -0.1039     0.6432 30.5962  -0.162 0.872692    
## rd.rpt                     1.7185     1.3072 35.9554   1.315 0.196956    
## grpd.ungr_lhd.rhd         -0.4688     1.4848 40.3254  -0.316 0.753815    
## grpd.ungr_cor.inc          9.5546     1.6008 26.7443   5.969 2.39e-06 ***
## cor.inc_lhd.rhd           -1.7708     1.1995 39.0258  -1.476 0.147889    
## grpd.ungr_rd.rpt           0.2958     2.1587 34.1389   0.137 0.891830    
## lhd.rhd_rd.rpt            -4.2403     1.7142 37.2433  -2.474 0.018056 *  
## cor.inc_rd.rpt            -0.7354     1.2808 29.9248  -0.574 0.570159    
## grpd.ungr_lhd.rhd_rd.rpt  -2.1463     2.0879 34.7450  -1.028 0.311067    
## cor.inc_lhd.rhd_rd.rpt     4.3161     2.4791 34.0037   1.741 0.090731 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
# correlations seem okay for Subj
# but not perfect for items
    'item     (Intercept)      6.750   2.598                                       
              lhd.rhd          4.960   2.227    -1.00             #not optimal                             
              cor.inc          2.834   1.683     0.88 -0.92                        
              cor.inc_lhd.rhd  8.860   2.977     0.95 -0.97  0.99 #not optimal'
## [1] "item     (Intercept)      6.750   2.598                                       \n              lhd.rhd          4.960   2.227    -1.00             #not optimal                             \n              cor.inc          2.834   1.683     0.88 -0.92                        \n              cor.inc_lhd.rhd  8.860   2.977     0.95 -0.97  0.99 #not optimal"
## check rePCA output =================================
summary(rePCA(m0_FL2))
## $Subj
## Importance of components:
##                          [,1]   [,2]   [,3]    [,4]    [,5]     [,6]      [,7]
## Standard deviation     1.4261 0.8921 0.8329 0.52991 0.28957 0.001824 2.546e-05
## Proportion of Variance 0.5231 0.2047 0.1784 0.07222 0.02157 0.000000 0.000e+00
## Cumulative Proportion  0.5231 0.7278 0.9062 0.97843 1.00000 1.000000 1.000e+00
## 
## $item
## Importance of components:
##                          [,1]    [,2]      [,3]      [,4]
## Standard deviation     0.9041 0.16070 0.0003666 6.829e-05
## Proportion of Variance 0.9694 0.03063 0.0000000 0.000e+00
## Cumulative Proportion  0.9694 1.00000 1.0000000 1.000e+00
# 5 of 7 components for Subj | 2 of 4 components for item contribute (variance)
# try model reduction of randEff ---------------------------------
## start with ZERO CORRELATION model =================================
summary(zcp_FL2 <- lmer(s8reln2 ~ 
                          1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                          grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                          cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                          lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                          grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                          (1+cor.inc+grpd.ungr+rd.rpt+
                             grpd.ungr_cor.inc+
                             cor.inc_rd.rpt+
                             grpd.ungr_rd.rpt||Subj)+
                          (1+lhd.rhd+cor.inc+
                             cor.inc_lhd.rhd||item),
                        REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## s8reln2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + cor.inc_rd.rpt + grpd.ungr_rd.rpt ||  
##     Subj) + (1 + lhd.rhd + cor.inc + cor.inc_lhd.rhd || item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8373.4   8502.4  -4161.7   8323.4     1264 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2302 -0.5851  0.0076  0.5721  4.0602 
## 
## Random effects:
##  Groups   Name              Variance  Std.Dev. 
##  Subj     (Intercept)       1.596e+01 3.994e+00
##  Subj.1   cor.inc           0.000e+00 0.000e+00
##  Subj.2   grpd.ungr         7.066e+00 2.658e+00
##  Subj.3   rd.rpt            1.214e+01 3.485e+00
##  Subj.4   grpd.ungr_cor.inc 6.061e+01 7.785e+00
##  Subj.5   cor.inc_rd.rpt    9.206e+00 3.034e+00
##  Subj.6   grpd.ungr_rd.rpt  1.050e+01 3.240e+00
##  item     (Intercept)       9.079e+00 3.013e+00
##  item.1   lhd.rhd           2.419e+00 1.555e+00
##  item.2   cor.inc           1.152e+00 1.073e+00
##  item.3   cor.inc_lhd.rhd   3.534e-09 5.945e-05
##  Residual                   2.858e+01 5.346e+00
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)               43.31334    0.97786  53.26175  44.294  < 2e-16 ***
## grpd.ungr                  4.25572    1.39282  31.88133   3.055  0.00452 ** 
## lhd.rhd                    3.91615    1.55288  35.95425   2.522  0.01624 *  
## cor.inc                   -0.01469    0.54154  54.67901  -0.027  0.97845    
## rd.rpt                     1.43481    1.48291  37.65768   0.968  0.33944    
## grpd.ungr_lhd.rhd         -1.10350    1.42861  35.65921  -0.772  0.44495    
## grpd.ungr_cor.inc         10.09133    1.72258  33.76543   5.858 1.35e-06 ***
## cor.inc_lhd.rhd           -1.56023    0.98699 443.50904  -1.581  0.11464    
## grpd.ungr_rd.rpt          -0.81302    2.65111  26.98838  -0.307  0.76145    
## lhd.rhd_rd.rpt            -4.54409    1.77172  36.27850  -2.565  0.01460 *  
## cor.inc_rd.rpt             0.01200    1.17695  27.03668   0.010  0.99194    
## grpd.ungr_lhd.rhd_rd.rpt  -0.67822    2.31076  28.99594  -0.294  0.77123    
## cor.inc_lhd.rhd_rd.rpt     4.29135    2.18515  31.85738   1.964  0.05832 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
    'Random effects:
     Groups   Name            Variance  Std.Dev. 
     Subj     (Intercept)       1.596e+01 3.9944841
     Subj.1   cor.inc           0.000e+00 0.0000000 #no variance
     Subj.2   grpd.ungr         7.065e+00 2.6579804
     Subj.3   rd.rpt            1.214e+01 3.4847284
     Subj.4   grpd.ungr_cor.inc 6.059e+01 7.7837605
     Subj.5   cor.inc_rd.rpt    9.206e+00 3.0340750
     Subj.6   grpd.ungr_rd.rpt  1.050e+01 3.2398369
     item     (Intercept)       9.081e+00 3.0134095
     item.1   lhd.rhd           2.418e+00 1.5549853
     item.2   cor.inc           1.152e+00 1.0732468
     item.3   cor.inc_lhd.rhd   3.940e-08 0.0001985 #no variance
     Residual                   2.858e+01 5.3458264
    Number of obs: 1289, groups:  Subj, 32; item, 24'
## [1] "Random effects:\n     Groups   Name            Variance  Std.Dev. \n     Subj     (Intercept)       1.596e+01 3.9944841\n     Subj.1   cor.inc           0.000e+00 0.0000000 #no variance\n     Subj.2   grpd.ungr         7.065e+00 2.6579804\n     Subj.3   rd.rpt            1.214e+01 3.4847284\n     Subj.4   grpd.ungr_cor.inc 6.059e+01 7.7837605\n     Subj.5   cor.inc_rd.rpt    9.206e+00 3.0340750\n     Subj.6   grpd.ungr_rd.rpt  1.050e+01 3.2398369\n     item     (Intercept)       9.081e+00 3.0134095\n     item.1   lhd.rhd           2.418e+00 1.5549853\n     item.2   cor.inc           1.152e+00 1.0732468\n     item.3   cor.inc_lhd.rhd   3.940e-08 0.0001985 #no variance\n     Residual                   2.858e+01 5.3458264\n    Number of obs: 1289, groups:  Subj, 32; item, 24"
## ANOVA model comparison - initial and zcp model =================================
anova(m0_FL2, zcp_FL2)
# m0_FL2 has signif. better fit - removing all correlations too harsh
## back to initial model =================================
# try exclusion of randEff slopes with regards to zcp + initial model output
### FIRST #############################
# try exclusion of "cor.inc" as slope per Subj
summary(m1_FL2 <- lmer(s8reln2 ~ 
                         1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                         grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                         cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                         lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                         grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                         (1+grpd.ungr+rd.rpt+
                            grpd.ungr_cor.inc+
                            cor.inc_rd.rpt+
                            grpd.ungr_rd.rpt|Subj)+
                         (1+lhd.rhd+cor.inc+
                            cor.inc_lhd.rhd|item),
                       REML = FALSE, pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## s8reln2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + grpd.ungr + rd.rpt + grpd.ungr_cor.inc +  
##     cor.inc_rd.rpt + grpd.ungr_rd.rpt | Subj) + (1 + lhd.rhd +  
##     cor.inc + cor.inc_lhd.rhd | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8352.4   8584.6  -4131.2   8262.4     1244 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2292 -0.6067 -0.0094  0.5770  4.0329 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)       16.212   4.026                                 
##           grpd.ungr          7.284   2.699    -0.73                        
##           rd.rpt            11.654   3.414    -0.01  0.05                  
##           grpd.ungr_cor.inc 59.585   7.719    -0.19  0.08 -0.54            
##           cor.inc_rd.rpt    10.642   3.262     0.12 -0.01 -0.09  0.20      
##           grpd.ungr_rd.rpt  10.048   3.170    -0.28 -0.10 -0.11  0.21 -0.87
##  item     (Intercept)        6.724   2.593                                 
##           lhd.rhd            4.784   2.187    -1.00                        
##           cor.inc            2.498   1.580     0.93 -0.93                  
##           cor.inc_lhd.rhd    8.017   2.831     0.98 -0.98  0.99            
##  Residual                   28.009   5.292                                 
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)              43.23713    0.92665 51.23790  46.659  < 2e-16 ***
## grpd.ungr                 4.10716    1.20773 36.79100   3.401  0.00163 ** 
## lhd.rhd                   3.84004    1.57169 37.67393   2.443  0.01936 *  
## cor.inc                   0.15422    0.55987 36.89006   0.275  0.78450    
## rd.rpt                    1.34602    1.32563 37.57666   1.015  0.31642    
## grpd.ungr_lhd.rhd        -0.75095    1.50785 41.50291  -0.498  0.62109    
## grpd.ungr_cor.inc        10.15091    1.69708 26.76707   5.981  2.3e-06 ***
## cor.inc_lhd.rhd          -1.66840    1.08237 52.68549  -1.541  0.12920    
## grpd.ungr_rd.rpt          0.05174    2.19425 36.05503   0.024  0.98132    
## lhd.rhd_rd.rpt           -4.50130    1.71621 37.79469  -2.623  0.01250 *  
## cor.inc_rd.rpt            0.38232    1.26105 29.98829   0.303  0.76385    
## grpd.ungr_lhd.rhd_rd.rpt -1.50242    2.09139 35.17177  -0.718  0.47726    
## cor.inc_lhd.rhd_rd.rpt    5.25565    2.44282 37.11923   2.151  0.03801 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
# Model failed to converge with 1 negative eigenvalue: -7.0e-01
### SECOND #############################
# back to initial model
# try exclusion of "cor.inc_lhd.rhd" as slope per item
summary(m2_FL2 <- lmer(s8reln2 ~ 
                         1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                         grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                         cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                         lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                         grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                         (1+cor.inc+grpd.ungr+rd.rpt+
                            grpd.ungr_cor.inc+
                            cor.inc_rd.rpt+
                            grpd.ungr_rd.rpt|Subj)+
                         (1+lhd.rhd+cor.inc|item),
                       REML = FALSE, pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## s8reln2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + cor.inc_rd.rpt + grpd.ungr_rd.rpt | Subj) +  
##     (1 + lhd.rhd + cor.inc | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8365.0   8612.7  -4134.5   8269.0     1241 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2583 -0.6091  0.0013  0.5985  3.9704 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)       17.448   4.177                                 
##           cor.inc            2.454   1.567    -0.43                        
##           grpd.ungr          6.652   2.579    -0.75  0.21                  
##           rd.rpt            10.095   3.177     0.02  0.47 -0.06            
##           grpd.ungr_cor.inc 49.346   7.025    -0.12 -0.84  0.23 -0.43      
##           cor.inc_rd.rpt    10.664   3.266     0.09 -0.07  0.18 -0.09  0.02
##           grpd.ungr_rd.rpt  10.598   3.256    -0.16 -0.25 -0.20 -0.23  0.35
##  item     (Intercept)        7.427   2.725                                 
##           lhd.rhd            2.137   1.462    -0.98                        
##           cor.inc            1.300   1.140     0.97 -0.92                  
##  Residual                   28.252   5.315                                 
##       
##       
##       
##       
##       
##       
##       
##  -0.88
##       
##       
##       
##       
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)               43.4684     0.9614 49.1543  45.215  < 2e-16 ***
## grpd.ungr                  4.5384     1.2737 31.0033   3.563  0.00121 ** 
## lhd.rhd                    3.8720     1.5907 32.2861   2.434  0.02063 *  
## cor.inc                   -0.1135     0.5848 32.5070  -0.194  0.84736    
## rd.rpt                     1.6312     1.3516 35.7998   1.207  0.23539    
## grpd.ungr_lhd.rhd         -0.8135     1.3748 35.5770  -0.592  0.55775    
## grpd.ungr_cor.inc          9.9653     1.5742 26.8695   6.330  9.1e-07 ***
## cor.inc_lhd.rhd           -1.5300     1.0086 54.7002  -1.517  0.13502    
## grpd.ungr_rd.rpt          -0.2401     2.3384 30.5457  -0.103  0.91889    
## lhd.rhd_rd.rpt            -4.3291     1.5665 32.0087  -2.764  0.00940 ** 
## cor.inc_rd.rpt            -0.4748     1.1800 30.6689  -0.402  0.69018    
## grpd.ungr_lhd.rhd_rd.rpt  -1.8997     2.1195 29.9436  -0.896  0.37725    
## cor.inc_lhd.rhd_rd.rpt     4.7796     2.1681 31.2683   2.204  0.03498 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
# Model failed to converge with 1 negative eigenvalue: -4.8e-03
### THIRD #############################
# back to initial model
# try exclude "lhd.rhd" as slope for item - had corr of -1 in "m0_FL2"
summary(m3_FL2 <- lmer(s8reln2 ~ 
                         1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                         grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                         cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                         lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                         grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                         (1+cor.inc+grpd.ungr+rd.rpt+
                            grpd.ungr_cor.inc+
                            cor.inc_rd.rpt+
                            grpd.ungr_rd.rpt|Subj)+
                         (1+cor.inc+
                            cor.inc_lhd.rhd|item),
                       REML = FALSE, pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## s8reln2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + cor.inc_rd.rpt + grpd.ungr_rd.rpt | Subj) +  
##     (1 + cor.inc + cor.inc_lhd.rhd | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8384.1   8631.9  -4144.1   8288.1     1241 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3323 -0.5920  0.0206  0.5999  4.0344 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)       17.3024  4.1596                                
##           cor.inc            2.5649  1.6015   -0.38                        
##           grpd.ungr          6.7412  2.5964   -0.76  0.14                  
##           rd.rpt             9.9679  3.1572    0.00  0.43 -0.05            
##           grpd.ungr_cor.inc 45.5727  6.7508   -0.12 -0.86  0.24 -0.39      
##           cor.inc_rd.rpt    10.4548  3.2334    0.15 -0.08  0.14 -0.07 -0.01
##           grpd.ungr_rd.rpt  10.6163  3.2583   -0.13 -0.28 -0.20 -0.25  0.38
##  item     (Intercept)        8.4238  2.9024                                
##           cor.inc            1.5891  1.2606   0.76                         
##           cor.inc_lhd.rhd    0.9492  0.9743   0.39  0.90                   
##  Residual                   28.6977  5.3570                                
##       
##       
##       
##       
##       
##       
##       
##  -0.88
##       
##       
##       
##       
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)               43.4922     0.9811 49.8430  44.330  < 2e-16 ***
## grpd.ungr                  4.4983     1.3402 30.4766   3.357  0.00213 ** 
## lhd.rhd                    3.9094     1.5588 30.2495   2.508  0.01773 *  
## cor.inc                   -0.1059     0.6006 28.2418  -0.176  0.86133    
## rd.rpt                     1.5621     1.4117 34.7016   1.107  0.27611    
## grpd.ungr_lhd.rhd         -0.8444     1.2487 29.7290  -0.676  0.50412    
## grpd.ungr_cor.inc         10.2051     1.5421 26.0626   6.618 5.04e-07 ***
## cor.inc_lhd.rhd           -1.3849     1.0358 44.8828  -1.337  0.18796    
## grpd.ungr_rd.rpt          -0.4127     2.5300 28.2477  -0.163  0.87158    
## lhd.rhd_rd.rpt            -4.4772     1.4662 26.0444  -3.054  0.00516 ** 
## cor.inc_rd.rpt            -0.5320     1.2029 27.4099  -0.442  0.66177    
## grpd.ungr_lhd.rhd_rd.rpt  -1.3452     1.8993 32.8101  -0.708  0.48379    
## cor.inc_lhd.rhd_rd.rpt     4.7692     2.2107 29.2855   2.157  0.03932 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular / Model failed to converge with 1 negative eigenvalue: -1.3e-03
# FINAL model ---------------------------------
## use initial model as final model =================================
# all model reduction steps lead to convergence fails
## initial model with factors and * for fixed effects =================================
# use REML=TRUE
summary(FL2 <- lmer(s8reln2 ~ 
                      condition*accuracy+condition*group*exp+
                      group*exp*accuracy+
                      (1+cor.inc+grpd.ungr+rd.rpt+
                         grpd.ungr_cor.inc+
                         cor.inc_rd.rpt+
                         grpd.ungr_rd.rpt|Subj)+
                      (1+lhd.rhd+cor.inc+
                         cor.inc_lhd.rhd|item),
                    REML = TRUE, pDat_mm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: s8reln2 ~ condition * accuracy + condition * group * exp + group *  
##     exp * accuracy + (1 + cor.inc + grpd.ungr + rd.rpt + grpd.ungr_cor.inc +  
##     cor.inc_rd.rpt + grpd.ungr_rd.rpt | Subj) + (1 + lhd.rhd +  
##     cor.inc + cor.inc_lhd.rhd | item)
##    Data: pDat_mm
## 
## REML criterion at convergence: 8222.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2224 -0.6160 -0.0031  0.5752  4.0520 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)       18.900   4.347                                 
##           cor.inc            2.964   1.722    -0.47                        
##           grpd.ungr          7.277   2.698    -0.72  0.29                  
##           rd.rpt            10.781   3.283     0.05  0.47 -0.07            
##           grpd.ungr_cor.inc 51.359   7.167    -0.13 -0.80  0.18 -0.43      
##           cor.inc_rd.rpt    13.132   3.624    -0.02 -0.02  0.18 -0.10  0.05
##           grpd.ungr_rd.rpt  12.222   3.496    -0.13 -0.31 -0.19 -0.24  0.39
##  item     (Intercept)        7.814   2.795                                 
##           lhd.rhd            5.639   2.375    -0.99                        
##           cor.inc            3.376   1.837     0.83 -0.88                  
##           cor.inc_lhd.rhd   10.222   3.197     0.92 -0.96  0.98            
##  Residual                   27.752   5.268                                 
##       
##       
##       
##       
##       
##       
##       
##  -0.84
##       
##       
##       
##       
##       
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                        Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)             43.4667     0.9946 45.0080  43.701  < 2e-16 ***
## condition1               4.7002     1.2599 31.2660   3.731  0.00076 ***
## accuracy1               -0.1078     0.6757 28.6669  -0.160  0.87434    
## group1                   3.9782     1.6939 32.7308   2.349  0.02504 *  
## exp1                     1.7214     1.3872 31.4757   1.241  0.22378    
## condition1:accuracy1     9.5091     1.6415 25.8839   5.793 4.28e-06 ***
## condition1:group1       -0.4451     1.5410 38.1461  -0.289  0.77426    
## condition1:exp1          0.1995     2.2935 28.0783   0.087  0.93132    
## group1:exp1             -4.2525     1.7861 34.6633  -2.381  0.02290 *  
## accuracy1:group1        -1.7256     1.2535 34.8840  -1.377  0.17738    
## accuracy1:exp1          -0.7481     1.3588 22.3199  -0.551  0.58740    
## condition1:group1:exp1  -2.2163     2.1720 30.4381  -1.020  0.31558    
## accuracy1:group1:exp1    4.3406     2.6122 24.0041   1.662  0.10959    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
save(FL2, file = "FL2.Rdata")
# *Nice table output* ---------------------------------
# sjPlot::tab_model((FL2), show.se = TRUE, show.stat = TRUE, collapse.ci=TRUE, file="FL2.doc")
# clean workspace ---------------------------------
# ls()
rm(m0_FL2, m1_FL2, m2_FL2, m3_FL2, zcp_FL2)

4.3 Model reduction - f0-range (F02)

# all "REML = FALSE" to use ANOVA model comparison later on
# INITIAL model with indicator variables ---------------------------------
summary(m0_F02 <- lmer(f0_range2 ~ 
                         1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                         grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                         cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                         lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                         grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                         (1+cor.inc+grpd.ungr+rd.rpt+
                            grpd.ungr_cor.inc+
                            cor.inc_rd.rpt+
                            grpd.ungr_rd.rpt|Subj)+
                         (1+lhd.rhd+cor.inc+
                            cor.inc_lhd.rhd|item),
                       REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## f0_range2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + cor.inc_rd.rpt + grpd.ungr_rd.rpt | Subj) +  
##     (1 + lhd.rhd + cor.inc + cor.inc_lhd.rhd | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   5955.8   6222.2  -2925.9   5851.8     1189 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4364 -0.4503 -0.0323  0.3961  6.8894 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)       3.05305  1.7473                                
##           cor.inc           0.06531  0.2556    1.00                        
##           grpd.ungr         1.35506  1.1641    0.48  0.48                  
##           rd.rpt            1.75924  1.3264   -0.12 -0.12  0.09            
##           grpd.ungr_cor.inc 7.10784  2.6661    0.49  0.49 -0.13 -0.37      
##           cor.inc_rd.rpt    0.93138  0.9651    0.83  0.83  0.55 -0.20  0.15
##           grpd.ungr_rd.rpt  2.31544  1.5217    0.01  0.01  0.60  0.81 -0.57
##  item     (Intercept)       0.03354  0.1831                                
##           lhd.rhd           0.52521  0.7247   -0.86                        
##           cor.inc           0.12813  0.3580    0.39 -0.81                  
##           cor.inc_lhd.rhd   0.40600  0.6372    0.92 -0.59  0.01            
##  Residual                   5.43790  2.3319                                
##       
##       
##       
##       
##       
##       
##       
##   0.09
##       
##       
##       
##       
##       
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                5.3222     0.3313 32.2559  16.064  < 2e-16 ***
## grpd.ungr                  2.0991     0.2890 28.6121   7.263 5.79e-08 ***
## lhd.rhd                    1.1229     0.6330 37.7712   1.774 0.084135 .  
## cor.inc                    0.3210     0.2230 77.4360   1.440 0.154028    
## rd.rpt                     2.3654     0.3295 29.4105   7.178 6.16e-08 ***
## grpd.ungr_lhd.rhd          0.2354     0.6014 35.0779   0.391 0.697878    
## grpd.ungr_cor.inc          3.0438     0.6349 26.6829   4.794 5.45e-05 ***
## cor.inc_lhd.rhd           -0.3729     0.4349 55.4545  -0.858 0.394835    
## grpd.ungr_rd.rpt           0.9006     0.4624 37.1558   1.948 0.059030 .  
## lhd.rhd_rd.rpt             2.9512     0.6819 36.1888   4.328 0.000114 ***
## cor.inc_rd.rpt             0.2083     0.4645 27.9049   0.448 0.657271    
## grpd.ungr_lhd.rhd_rd.rpt   1.8301     0.9141 39.4978   2.002 0.052184 .  
## cor.inc_lhd.rhd_rd.rpt     1.2205     0.9171 23.6050   1.331 0.195959    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
# Model failed to converge with 2 negative eigenvalues: -1.5e-03 -1.6e+01
summary(rePCA(m0_F02))
## $Subj
## Importance of components:
##                          [,1]   [,2]   [,3]   [,4]    [,5]      [,6]      [,7]
## Standard deviation     1.3133 0.8990 0.6007 0.3659 0.15069 4.165e-05 1.417e-17
## Proportion of Variance 0.5654 0.2650 0.1183 0.0439 0.00744 0.000e+00 0.000e+00
## Cumulative Proportion  0.5654 0.8304 0.9487 0.9926 1.00000 1.000e+00 1.000e+00
## 
## $item
## Importance of components:
##                          [,1]   [,2]      [,3] [,4]
## Standard deviation     0.3884 0.2238 0.0002056    0
## Proportion of Variance 0.7508 0.2492 0.0000000    0
## Cumulative Proportion  0.7508 1.0000 1.0000000    1
# 5 of 7 components for Subj | 2 of 4 components for item contribute (variance)
# try model reduction of randEff ---------------------------------
## start with ZERO CORRELATION model =================================
summary(zcp_F02 <- lmer(f0_range2 ~ 
                          1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                          grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                          cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                          lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                          grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                          (1+cor.inc+grpd.ungr+rd.rpt+
                             grpd.ungr_cor.inc+
                             cor.inc_rd.rpt+
                             grpd.ungr_rd.rpt||Subj)+
                          (1+lhd.rhd+cor.inc+
                             cor.inc_lhd.rhd||item),
                        REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## f0_range2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + cor.inc_rd.rpt + grpd.ungr_rd.rpt ||  
##     Subj) + (1 + lhd.rhd + cor.inc + cor.inc_lhd.rhd || item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   5949.1   6077.2  -2949.6   5899.1     1216 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0695 -0.4312 -0.0440  0.3984  6.7809 
## 
## Random effects:
##  Groups   Name              Variance  Std.Dev. 
##  Subj     (Intercept)       3.238e+00 1.799e+00
##  Subj.1   cor.inc           0.000e+00 0.000e+00
##  Subj.2   grpd.ungr         1.220e+00 1.104e+00
##  Subj.3   rd.rpt            1.469e+00 1.212e+00
##  Subj.4   grpd.ungr_cor.inc 5.502e+00 2.346e+00
##  Subj.5   cor.inc_rd.rpt    1.047e+00 1.023e+00
##  Subj.6   grpd.ungr_rd.rpt  4.089e-01 6.394e-01
##  item     (Intercept)       4.275e-02 2.067e-01
##  item.1   lhd.rhd           4.676e-01 6.838e-01
##  item.2   cor.inc           0.000e+00 0.000e+00
##  item.3   cor.inc_lhd.rhd   6.944e-09 8.333e-05
##  Residual                   5.631e+00 2.373e+00
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                5.3717     0.3417  35.2945  15.722  < 2e-16 ***
## grpd.ungr                  2.0899     0.2842  31.3584   7.354 2.61e-08 ***
## lhd.rhd                    1.6858     0.6924  36.9313   2.435 0.019847 *  
## cor.inc                    0.4174     0.2115 411.0690   1.974 0.049099 *  
## rd.rpt                     2.2444     0.3228  36.6579   6.954 3.43e-08 ***
## grpd.ungr_lhd.rhd          0.1203     0.5973  33.5177   0.201 0.841639    
## grpd.ungr_cor.inc          3.2277     0.5741  32.5001   5.622 3.09e-06 ***
## cor.inc_lhd.rhd           -0.2549     0.4215 430.0617  -0.605 0.545551    
## grpd.ungr_rd.rpt           0.7294     0.3836  23.4116   1.902 0.069589 .  
## lhd.rhd_rd.rpt             2.5498     0.6815  39.9618   3.741 0.000575 ***
## cor.inc_rd.rpt             0.2252     0.4422  29.2211   0.509 0.614403    
## grpd.ungr_lhd.rhd_rd.rpt   1.3568     0.8804  23.7869   1.541 0.136478    
## cor.inc_lhd.rhd_rd.rpt     1.2258     0.8847  29.2133   1.386 0.176386    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
summary(zcp_F02)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## f0_range2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + cor.inc_rd.rpt + grpd.ungr_rd.rpt ||  
##     Subj) + (1 + lhd.rhd + cor.inc + cor.inc_lhd.rhd || item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   5949.1   6077.2  -2949.6   5899.1     1216 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0695 -0.4312 -0.0440  0.3984  6.7809 
## 
## Random effects:
##  Groups   Name              Variance  Std.Dev. 
##  Subj     (Intercept)       3.238e+00 1.799e+00
##  Subj.1   cor.inc           0.000e+00 0.000e+00
##  Subj.2   grpd.ungr         1.220e+00 1.104e+00
##  Subj.3   rd.rpt            1.469e+00 1.212e+00
##  Subj.4   grpd.ungr_cor.inc 5.502e+00 2.346e+00
##  Subj.5   cor.inc_rd.rpt    1.047e+00 1.023e+00
##  Subj.6   grpd.ungr_rd.rpt  4.089e-01 6.394e-01
##  item     (Intercept)       4.275e-02 2.067e-01
##  item.1   lhd.rhd           4.676e-01 6.838e-01
##  item.2   cor.inc           0.000e+00 0.000e+00
##  item.3   cor.inc_lhd.rhd   6.944e-09 8.333e-05
##  Residual                   5.631e+00 2.373e+00
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                5.3717     0.3417  35.2945  15.722  < 2e-16 ***
## grpd.ungr                  2.0899     0.2842  31.3584   7.354 2.61e-08 ***
## lhd.rhd                    1.6858     0.6924  36.9313   2.435 0.019847 *  
## cor.inc                    0.4174     0.2115 411.0690   1.974 0.049099 *  
## rd.rpt                     2.2444     0.3228  36.6579   6.954 3.43e-08 ***
## grpd.ungr_lhd.rhd          0.1203     0.5973  33.5177   0.201 0.841639    
## grpd.ungr_cor.inc          3.2277     0.5741  32.5001   5.622 3.09e-06 ***
## cor.inc_lhd.rhd           -0.2549     0.4215 430.0617  -0.605 0.545551    
## grpd.ungr_rd.rpt           0.7294     0.3836  23.4116   1.902 0.069589 .  
## lhd.rhd_rd.rpt             2.5498     0.6815  39.9618   3.741 0.000575 ***
## cor.inc_rd.rpt             0.2252     0.4422  29.2211   0.509 0.614403    
## grpd.ungr_lhd.rhd_rd.rpt   1.3568     0.8804  23.7869   1.541 0.136478    
## cor.inc_lhd.rhd_rd.rpt     1.2258     0.8847  29.2133   1.386 0.176386    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
    'Random effects:
     Groups   Name            Variance  Std.Dev. 
     Subj     (Intercept)       3.238e+00 1.7993344
     Subj.1   cor.inc           0.000e+00 0.0000000 # no variance
     Subj.2   grpd.ungr         1.220e+00 1.1043709
     Subj.3   rd.rpt            1.469e+00 1.2118683
     Subj.4   grpd.ungr_cor.inc 5.502e+00 2.3457306
     Subj.5   cor.inc_rd.rpt    1.047e+00 1.0231915 # mean var. + sd almost equal
     Subj.6   grpd.ungr_rd.rpt  4.087e-01 0.6392808 # mean var. smaller than sd
     item     (Intercept)       4.275e-02 0.2067492
     item.1   lhd.rhd           4.676e-01 0.6838274
     item.2   cor.inc           8.172e-08 0.0002859 # no variance
     item.3   cor.inc_lhd.rhd   0.000e+00 0.0000000 # no variance
     Residual                   5.631e+00 2.3728895
    Number of obs: 1241, groups:  Subj, 32; item, 24'
## [1] "Random effects:\n     Groups   Name            Variance  Std.Dev. \n     Subj     (Intercept)       3.238e+00 1.7993344\n     Subj.1   cor.inc           0.000e+00 0.0000000 # no variance\n     Subj.2   grpd.ungr         1.220e+00 1.1043709\n     Subj.3   rd.rpt            1.469e+00 1.2118683\n     Subj.4   grpd.ungr_cor.inc 5.502e+00 2.3457306\n     Subj.5   cor.inc_rd.rpt    1.047e+00 1.0231915 # mean var. + sd almost equal\n     Subj.6   grpd.ungr_rd.rpt  4.087e-01 0.6392808 # mean var. smaller than sd\n     item     (Intercept)       4.275e-02 0.2067492\n     item.1   lhd.rhd           4.676e-01 0.6838274\n     item.2   cor.inc           8.172e-08 0.0002859 # no variance\n     item.3   cor.inc_lhd.rhd   0.000e+00 0.0000000 # no variance\n     Residual                   5.631e+00 2.3728895\n    Number of obs: 1241, groups:  Subj, 32; item, 24"
## ANOVA model comparison - initial - zcp model =================================
anova(m0_F02, zcp_F02)
# initial model m0_F02 significantly better variance resolution
# BUT zcp_F02 has smaller AIC - therefore better model fit
# reasonable to look further into zcp model
## iterative reduction of randEff =================================
### FIRST #############################
# try exclusion of "cor.inc" as slope per Subj
summary(zcp1_F02 <- lmer(f0_range2 ~ 
                           1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                           grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                           cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                           lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                           grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                           (1+grpd.ungr+rd.rpt+
                              grpd.ungr_cor.inc+
                              cor.inc_rd.rpt+
                              grpd.ungr_rd.rpt||Subj)+
                           (1+lhd.rhd+cor.inc+
                              cor.inc_lhd.rhd||item),
                         REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## f0_range2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + grpd.ungr + rd.rpt + grpd.ungr_cor.inc +  
##     cor.inc_rd.rpt + grpd.ungr_rd.rpt || Subj) + (1 + lhd.rhd +  
##     cor.inc + cor.inc_lhd.rhd || item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   5947.1   6070.1  -2949.6   5899.1     1217 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0695 -0.4312 -0.0440  0.3985  6.7809 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev.
##  Subj     (Intercept)       3.23785  1.7994  
##  Subj.1   grpd.ungr         1.21953  1.1043  
##  Subj.2   rd.rpt            1.46871  1.2119  
##  Subj.3   grpd.ungr_cor.inc 5.50201  2.3456  
##  Subj.4   cor.inc_rd.rpt    1.04687  1.0232  
##  Subj.5   grpd.ungr_rd.rpt  0.40876  0.6393  
##  item     (Intercept)       0.04275  0.2068  
##  item.1   lhd.rhd           0.46757  0.6838  
##  item.2   cor.inc           0.00000  0.0000  
##  item.3   cor.inc_lhd.rhd   0.00000  0.0000  
##  Residual                   5.63060  2.3729  
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                5.3717     0.3417  35.2937  15.722  < 2e-16 ***
## grpd.ungr                  2.0899     0.2842  31.3596   7.354 2.61e-08 ***
## lhd.rhd                    1.6858     0.6924  36.9304   2.435 0.019848 *  
## cor.inc                    0.4174     0.2115 411.0560   1.974 0.049099 *  
## rd.rpt                     2.2444     0.3228  36.6586   6.954 3.43e-08 ***
## grpd.ungr_lhd.rhd          0.1203     0.5973  33.5181   0.201 0.841642    
## grpd.ungr_cor.inc          3.2277     0.5741  32.5008   5.622 3.09e-06 ***
## cor.inc_lhd.rhd           -0.2550     0.4214 430.0484  -0.605 0.545537    
## grpd.ungr_rd.rpt           0.7294     0.3836  23.4103   1.902 0.069589 .  
## lhd.rhd_rd.rpt             2.5498     0.6815  39.9620   3.741 0.000575 ***
## cor.inc_rd.rpt             0.2252     0.4422  29.2201   0.509 0.614390    
## grpd.ungr_lhd.rhd_rd.rpt   1.3569     0.8804  23.7853   1.541 0.136469    
## cor.inc_lhd.rhd_rd.rpt     1.2258     0.8847  29.2124   1.386 0.176374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
    'Random effects:
     Groups   Name            Variance Std.Dev.
     Subj     (Intercept)       3.23828  1.7995  
     Subj.1   grpd.ungr         1.21951  1.1043  
     Subj.2   rd.rpt            1.46867  1.2119  
     Subj.3   grpd.ungr_cor.inc 5.50127  2.3455  
     Subj.4   cor.inc_rd.rpt    1.04691  1.0232 # mean var. + sd almost equal 
     Subj.5   grpd.ungr_rd.rpt  0.40879  0.6394 # mean var. smaller than sd 
     item     (Intercept)       0.04276  0.2068  
     item.1   lhd.rhd           0.46757  0.6838  
     item.2   cor.inc           0.00000  0.0000 # no variance 
     item.3   cor.inc_lhd.rhd   0.00000  0.0000 # no variance
     Residual                   5.63059  2.3729  
    Number of obs: 1241, groups:  Subj, 32; item, 24'
## [1] "Random effects:\n     Groups   Name            Variance Std.Dev.\n     Subj     (Intercept)       3.23828  1.7995  \n     Subj.1   grpd.ungr         1.21951  1.1043  \n     Subj.2   rd.rpt            1.46867  1.2119  \n     Subj.3   grpd.ungr_cor.inc 5.50127  2.3455  \n     Subj.4   cor.inc_rd.rpt    1.04691  1.0232 # mean var. + sd almost equal \n     Subj.5   grpd.ungr_rd.rpt  0.40879  0.6394 # mean var. smaller than sd \n     item     (Intercept)       0.04276  0.2068  \n     item.1   lhd.rhd           0.46757  0.6838  \n     item.2   cor.inc           0.00000  0.0000 # no variance \n     item.3   cor.inc_lhd.rhd   0.00000  0.0000 # no variance\n     Residual                   5.63059  2.3729  \n    Number of obs: 1241, groups:  Subj, 32; item, 24"
summary(rePCA(zcp1_F02))
## $Subj
## Importance of components:
##                          [,1]   [,2]   [,3]    [,4]    [,5]    [,6]
## Standard deviation     0.9885 0.7583 0.5107 0.46539 0.43119 0.26944
## Proportion of Variance 0.4270 0.2513 0.1140 0.09466 0.08125 0.03173
## Cumulative Proportion  0.4270 0.6784 0.7924 0.88702 0.96827 1.00000
## 
## $item
## Importance of components:
##                          [,1]    [,2] [,3] [,4]
## Standard deviation     0.2882 0.08714    0    0
## Proportion of Variance 0.9162 0.08378    0    0
## Cumulative Proportion  0.9162 1.00000    1    1
# components for Subj seem okay | 2 of 4 components for item have zero variance
### ANOVA model comparison - reduced zcp - zcp model #############################
anova(zcp_F02, zcp1_F02)
# not significant, AIC smaller for zcp1_F02 - preferred
### SECOND #############################
# try also exclusion also of "cor.inc_lhd.rhd" as slope per item
summary(zcp2_F02 <- lmer(f0_range2 ~ 
                           1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                           grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                           cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                           lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                           grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                           (1+grpd.ungr+rd.rpt+
                              grpd.ungr_cor.inc+
                              cor.inc_rd.rpt+
                              grpd.ungr_rd.rpt||Subj)+
                           (1+lhd.rhd+cor.inc||item),
                         REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## f0_range2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + grpd.ungr + rd.rpt + grpd.ungr_cor.inc +  
##     cor.inc_rd.rpt + grpd.ungr_rd.rpt || Subj) + (1 + lhd.rhd +  
##     cor.inc || item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   5945.1   6062.9  -2949.6   5899.1     1218 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0695 -0.4312 -0.0440  0.3984  6.7809 
## 
## Random effects:
##  Groups   Name              Variance  Std.Dev. 
##  Subj     (Intercept)       3.238e+00 1.799e+00
##  Subj.1   grpd.ungr         1.219e+00 1.104e+00
##  Subj.2   rd.rpt            1.469e+00 1.212e+00
##  Subj.3   grpd.ungr_cor.inc 5.503e+00 2.346e+00
##  Subj.4   cor.inc_rd.rpt    1.047e+00 1.023e+00
##  Subj.5   grpd.ungr_rd.rpt  4.090e-01 6.395e-01
##  item     (Intercept)       4.274e-02 2.067e-01
##  item.1   lhd.rhd           4.675e-01 6.838e-01
##  item.2   cor.inc           4.418e-09 6.647e-05
##  Residual                   5.631e+00 2.373e+00
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                5.3717     0.3417  35.2942  15.722  < 2e-16 ***
## grpd.ungr                  2.0899     0.2842  31.3593   7.354 2.61e-08 ***
## lhd.rhd                    1.6858     0.6924  36.9310   2.435 0.019847 *  
## cor.inc                    0.4174     0.2115 411.0646   1.974 0.049091 *  
## rd.rpt                     2.2444     0.3228  36.6570   6.953 3.43e-08 ***
## grpd.ungr_lhd.rhd          0.1203     0.5973  33.5188   0.201 0.841652    
## grpd.ungr_cor.inc          3.2277     0.5741  32.4985   5.622 3.09e-06 ***
## cor.inc_lhd.rhd           -0.2549     0.4215 430.0566  -0.605 0.545597    
## grpd.ungr_rd.rpt           0.7294     0.3836  23.4124   1.902 0.069593 .  
## lhd.rhd_rd.rpt             2.5498     0.6815  39.9611   3.741 0.000575 ***
## cor.inc_rd.rpt             0.2252     0.4422  29.2213   0.509 0.614389    
## grpd.ungr_lhd.rhd_rd.rpt   1.3568     0.8804  23.7880   1.541 0.136493    
## cor.inc_lhd.rhd_rd.rpt     1.2258     0.8847  29.2135   1.386 0.176380    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
    'Random effects:
     Groups   Name            Variance Std.Dev.
     Subj     (Intercept)       3.23777  1.7994  
     Subj.1   grpd.ungr         1.21950  1.1043  
     Subj.2   rd.rpt            1.46869  1.2119  
     Subj.3   grpd.ungr_cor.inc 5.50204  2.3456  
     Subj.4   cor.inc_rd.rpt    1.04707  1.0233 # mean var. + sd almost equal  
     Subj.5   grpd.ungr_rd.rpt  0.40879  0.6394 # mean var. smaller than sd
     item     (Intercept)       0.04275  0.2068  
     item.1   lhd.rhd           0.46755  0.6838  
     item.2   cor.inc           0.00000  0.0000 # no variance
     Residual                   5.63061  2.3729  
    Number of obs: 1241, groups:  Subj, 32; item, 24'
## [1] "Random effects:\n     Groups   Name            Variance Std.Dev.\n     Subj     (Intercept)       3.23777  1.7994  \n     Subj.1   grpd.ungr         1.21950  1.1043  \n     Subj.2   rd.rpt            1.46869  1.2119  \n     Subj.3   grpd.ungr_cor.inc 5.50204  2.3456  \n     Subj.4   cor.inc_rd.rpt    1.04707  1.0233 # mean var. + sd almost equal  \n     Subj.5   grpd.ungr_rd.rpt  0.40879  0.6394 # mean var. smaller than sd\n     item     (Intercept)       0.04275  0.2068  \n     item.1   lhd.rhd           0.46755  0.6838  \n     item.2   cor.inc           0.00000  0.0000 # no variance\n     Residual                   5.63061  2.3729  \n    Number of obs: 1241, groups:  Subj, 32; item, 24"
summary(rePCA(zcp2_F02))
## $Subj
## Importance of components:
##                          [,1]   [,2]   [,3]    [,4]    [,5]    [,6]
## Standard deviation     0.9886 0.7583 0.5107 0.46538 0.43127 0.26951
## Proportion of Variance 0.4271 0.2513 0.1140 0.09464 0.08128 0.03174
## Cumulative Proportion  0.4271 0.6784 0.7923 0.88698 0.96826 1.00000
## 
## $item
## Importance of components:
##                          [,1]    [,2]      [,3]
## Standard deviation     0.2882 0.08713 2.801e-05
## Proportion of Variance 0.9162 0.08376 0.000e+00
## Cumulative Proportion  0.9162 1.00000 1.000e+00
# components for Subj seem okay - 1 of 3 components for item has zero variance
### ANOVA model comparison - reduced 2nd zcp - reduced 1st zcp model #############################
anova(zcp2_F02, zcp1_F02)
# not significant and AIC smaller for zcp2_F02 - preferred
### THIRD #############################
# try exclusion also of "cor.inc" as slope per item
summary(zcp3_F02 <- lmer(f0_range2 ~ 
                           1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                           grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                           cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                           lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                           grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                           (1+grpd.ungr+rd.rpt+
                              grpd.ungr_cor.inc+
                              cor.inc_rd.rpt+
                              grpd.ungr_rd.rpt||Subj)+
                           (1+lhd.rhd||item),
                         REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## f0_range2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + grpd.ungr + rd.rpt + grpd.ungr_cor.inc +  
##     cor.inc_rd.rpt + grpd.ungr_rd.rpt || Subj) + (1 + lhd.rhd ||      item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   5943.1   6055.8  -2949.6   5899.1     1219 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0695 -0.4312 -0.0440  0.3984  6.7809 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev.
##  Subj     (Intercept)       3.23798  1.7994  
##  Subj.1   grpd.ungr         1.21956  1.1043  
##  Subj.2   rd.rpt            1.46878  1.2119  
##  Subj.3   grpd.ungr_cor.inc 5.50252  2.3457  
##  Subj.4   cor.inc_rd.rpt    1.04718  1.0233  
##  Subj.5   grpd.ungr_rd.rpt  0.40874  0.6393  
##  item     (Intercept)       0.04274  0.2067  
##  item.1   lhd.rhd           0.46754  0.6838  
##  Residual                   5.63059  2.3729  
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                5.3717     0.3417  35.2911  15.721  < 2e-16 ***
## grpd.ungr                  2.0899     0.2842  31.3580   7.354 2.61e-08 ***
## lhd.rhd                    1.6858     0.6924  36.9278   2.435 0.019851 *  
## cor.inc                    0.4174     0.2115 411.0610   1.974 0.049095 *  
## rd.rpt                     2.2444     0.3228  36.6557   6.953 3.43e-08 ***
## grpd.ungr_lhd.rhd          0.1203     0.5973  33.5178   0.201 0.841660    
## grpd.ungr_cor.inc          3.2277     0.5741  32.4987   5.622 3.09e-06 ***
## cor.inc_lhd.rhd           -0.2549     0.4215 430.0530  -0.605 0.545577    
## grpd.ungr_rd.rpt           0.7294     0.3835  23.4097   1.902 0.069579 .  
## lhd.rhd_rd.rpt             2.5498     0.6815  39.9601   3.741 0.000575 ***
## cor.inc_rd.rpt             0.2252     0.4422  29.2210   0.509 0.614373    
## grpd.ungr_lhd.rhd_rd.rpt   1.3568     0.8804  23.7859   1.541 0.136476    
## cor.inc_lhd.rhd_rd.rpt     1.2258     0.8847  29.2132   1.386 0.176381    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    'Random effects:
     Groups   Name            Variance Std.Dev.
     Subj     (Intercept)       3.23788  1.7994  
     Subj.1   grpd.ungr         1.21961  1.1044  
     Subj.2   rd.rpt            1.46866  1.2119  
     Subj.3   grpd.ungr_cor.inc 5.50221  2.3457  
     Subj.4   cor.inc_rd.rpt    1.04710  1.0233 # mean var. + sd almost equal 
     Subj.5   grpd.ungr_rd.rpt  0.40869  0.6393 # mean var. smaller than sd 
     item     (Intercept)       0.04274  0.2067  
     item.1   lhd.rhd           0.46756  0.6838  
     Residual                   5.63060  2.3729  
    Number of obs: 1241, groups:  Subj, 32; item, 24'
## [1] "Random effects:\n     Groups   Name            Variance Std.Dev.\n     Subj     (Intercept)       3.23788  1.7994  \n     Subj.1   grpd.ungr         1.21961  1.1044  \n     Subj.2   rd.rpt            1.46866  1.2119  \n     Subj.3   grpd.ungr_cor.inc 5.50221  2.3457  \n     Subj.4   cor.inc_rd.rpt    1.04710  1.0233 # mean var. + sd almost equal \n     Subj.5   grpd.ungr_rd.rpt  0.40869  0.6393 # mean var. smaller than sd \n     item     (Intercept)       0.04274  0.2067  \n     item.1   lhd.rhd           0.46756  0.6838  \n     Residual                   5.63060  2.3729  \n    Number of obs: 1241, groups:  Subj, 32; item, 24"
summary(rePCA(zcp3_F02))
## $Subj
## Importance of components:
##                          [,1]   [,2]   [,3]    [,4]    [,5]    [,6]
## Standard deviation     0.9886 0.7583 0.5107 0.46540 0.43125 0.26943
## Proportion of Variance 0.4271 0.2513 0.1140 0.09465 0.08127 0.03172
## Cumulative Proportion  0.4271 0.6784 0.7923 0.88701 0.96828 1.00000
## 
## $item
## Importance of components:
##                          [,1]    [,2]
## Standard deviation     0.2882 0.08712
## Proportion of Variance 0.9162 0.08376
## Cumulative Proportion  0.9162 1.00000
# all components seem okay
### ANOVA model comparison - reduced 3rd zcp - reduced 2nd zcp model #############################
anova(zcp3_F02, zcp2_F02)
# not significant and AIC smaller for zcp3_F02 - preferred
### FOURTH #############################
# try exclusion also of "grpd.ungr_rd.rpt" as slope per Subj
summary(zcp4_F02 <- lmer(f0_range2 ~ 
                           1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                           grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                           cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                           lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                           grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                           (1+grpd.ungr+rd.rpt+
                              grpd.ungr_cor.inc+
                              cor.inc_rd.rpt||Subj)+
                           (1+lhd.rhd||item),
                         REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## f0_range2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + grpd.ungr + rd.rpt + grpd.ungr_cor.inc +  
##     cor.inc_rd.rpt || Subj) + (1 + lhd.rhd || item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   5941.5   6049.0  -2949.7   5899.5     1220 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0130 -0.4368 -0.0452  0.4072  6.7353 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev.
##  Subj     (Intercept)       3.2344   1.7984  
##  Subj.1   grpd.ungr         1.2075   1.0989  
##  Subj.2   rd.rpt            1.4482   1.2034  
##  Subj.3   grpd.ungr_cor.inc 5.4742   2.3397  
##  Subj.4   cor.inc_rd.rpt    1.1115   1.0543  
##  item     (Intercept)       0.0418   0.2045  
##  item.1   lhd.rhd           0.4660   0.6827  
##  Residual                   5.6543   2.3779  
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                5.3698     0.3414  35.3045  15.727  < 2e-16 ***
## grpd.ungr                  2.0785     0.2824  31.2622   7.359 2.63e-08 ***
## lhd.rhd                    1.6888     0.6920  36.9575   2.441 0.019573 *  
## cor.inc                    0.4239     0.2110 416.2915   2.009 0.045174 *  
## rd.rpt                     2.2337     0.3214  36.6877   6.949 3.46e-08 ***
## grpd.ungr_lhd.rhd          0.1079     0.5942  33.4156   0.182 0.857023    
## grpd.ungr_cor.inc          3.2464     0.5716  32.5122   5.680 2.61e-06 ***
## cor.inc_lhd.rhd           -0.2559     0.4206 433.2385  -0.608 0.543277    
## grpd.ungr_rd.rpt           0.7355     0.3617  37.1966   2.034 0.049169 *  
## lhd.rhd_rd.rpt             2.5438     0.6793  40.0039   3.745 0.000569 ***
## cor.inc_rd.rpt             0.2514     0.4421  28.8661   0.569 0.573949    
## grpd.ungr_lhd.rhd_rd.rpt   1.3749     0.8438  27.7374   1.629 0.114535    
## cor.inc_lhd.rhd_rd.rpt     1.2140     0.8844  28.8625   1.373 0.180424    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00466327 (tol = 0.002, component 1)
    'Random effects:
     Groups   Name            Variance Std.Dev.
     Subj     (Intercept)       3.23610  1.7989  
     Subj.1   grpd.ungr         1.20739  1.0988  
     Subj.2   rd.rpt            1.44618  1.2026  
     Subj.3   grpd.ungr_cor.inc 5.47603  2.3401  
     Subj.4   cor.inc_rd.rpt    1.11054  1.0538 # mean var. + sd almost equal 
     item     (Intercept)       0.04183  0.2045  
     item.1   lhd.rhd           0.46619  0.6828  
     Residual                   5.65429  2.3779  
    Number of obs: 1241, groups:  Subj, 32; item, 24'
## [1] "Random effects:\n     Groups   Name            Variance Std.Dev.\n     Subj     (Intercept)       3.23610  1.7989  \n     Subj.1   grpd.ungr         1.20739  1.0988  \n     Subj.2   rd.rpt            1.44618  1.2026  \n     Subj.3   grpd.ungr_cor.inc 5.47603  2.3401  \n     Subj.4   cor.inc_rd.rpt    1.11054  1.0538 # mean var. + sd almost equal \n     item     (Intercept)       0.04183  0.2045  \n     item.1   lhd.rhd           0.46619  0.6828  \n     Residual                   5.65429  2.3779  \n    Number of obs: 1241, groups:  Subj, 32; item, 24"
summary(rePCA(zcp4_F02))
## $Subj
## Importance of components:
##                          [,1]   [,2]   [,3]    [,4]    [,5]
## Standard deviation     0.9839 0.7563 0.5061 0.46212 0.44337
## Proportion of Variance 0.4388 0.2592 0.1161 0.09679 0.08909
## Cumulative Proportion  0.4388 0.6980 0.8141 0.91091 1.00000
## 
## $item
## Importance of components:
##                          [,1]    [,2]
## Standard deviation     0.2871 0.08598
## Proportion of Variance 0.9177 0.08231
## Cumulative Proportion  0.9177 1.00000
# all components seem okay
### ANOVA model comparison - reduced 4th zcp - reduced 3rd zcp model #############################
anova(zcp4_F02, zcp3_F02)
# not significant and AIC smaller for zcp4_F02 - preferred
### FIFTH #############################
# try exclusion also of "cor.inc_rd.rpt" as slope per Subj
summary(zcp5_F02 <- lmer(f0_range2 ~ 
                           1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                           grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                           cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                           lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                           grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                           (1+grpd.ungr+rd.rpt+
                              grpd.ungr_cor.inc||Subj)+
                           (1+lhd.rhd||item),
                         REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## f0_range2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + grpd.ungr + rd.rpt + grpd.ungr_cor.inc ||  
##     Subj) + (1 + lhd.rhd || item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   5940.3   6042.8  -2950.2   5900.3     1221 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0144 -0.4412 -0.0484  0.4113  6.7447 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev.
##  Subj     (Intercept)       3.29998  1.8166  
##  Subj.1   grpd.ungr         1.16990  1.0816  
##  Subj.2   rd.rpt            1.59736  1.2639  
##  Subj.3   grpd.ungr_cor.inc 5.49785  2.3447  
##  item     (Intercept)       0.04001  0.2000  
##  item.1   lhd.rhd           0.45381  0.6737  
##  Residual                   5.69228  2.3858  
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                 5.36753    0.34430   35.47415  15.590  < 2e-16 ***
## grpd.ungr                   2.06464    0.27855   31.26358   7.412 2.27e-08 ***
## lhd.rhd                     1.68259    0.69746   37.06045   2.412 0.020910 *  
## cor.inc                     0.40069    0.20595  925.92923   1.946 0.052006 .  
## rd.rpt                      2.28622    0.32548   36.38732   7.024 2.88e-08 ***
## grpd.ungr_lhd.rhd           0.09631    0.58684   33.40822   0.164 0.870634    
## grpd.ungr_cor.inc           3.24632    0.57027   32.65268   5.693 2.47e-06 ***
## cor.inc_lhd.rhd            -0.28707    0.41072  963.49607  -0.699 0.484764    
## grpd.ungr_rd.rpt            0.75711    0.35575   36.10486   2.128 0.040213 *  
## lhd.rhd_rd.rpt              2.59721    0.68657   39.77712   3.783 0.000511 ***
## cor.inc_rd.rpt              0.21811    0.39229 1126.74603   0.556 0.578332    
## grpd.ungr_lhd.rhd_rd.rpt    1.51546    0.83233   27.13287   1.821 0.079694 .  
## cor.inc_lhd.rhd_rd.rpt      1.23852    0.78475 1121.37900   1.578 0.114794    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    'Random effects:
     Groups   Name            Variance Std.Dev.
     Subj     (Intercept)       3.29999  1.8166  
     Subj.1   grpd.ungr         1.16990  1.0816  
     Subj.2   rd.rpt            1.59736  1.2639  
     Subj.3   grpd.ungr_cor.inc 5.49785  2.3447  
     item     (Intercept)       0.04001  0.2000  
     item.1   lhd.rhd           0.45381  0.6737  
     Residual                   5.69227  2.3858  
    Number of obs: 1241, groups:  Subj, 32; item, 24'
## [1] "Random effects:\n     Groups   Name            Variance Std.Dev.\n     Subj     (Intercept)       3.29999  1.8166  \n     Subj.1   grpd.ungr         1.16990  1.0816  \n     Subj.2   rd.rpt            1.59736  1.2639  \n     Subj.3   grpd.ungr_cor.inc 5.49785  2.3447  \n     item     (Intercept)       0.04001  0.2000  \n     item.1   lhd.rhd           0.45381  0.6737  \n     Residual                   5.69227  2.3858  \n    Number of obs: 1241, groups:  Subj, 32; item, 24"
summary(rePCA(zcp5_F02))
## $Subj
## Importance of components:
##                          [,1]   [,2]   [,3]   [,4]
## Standard deviation     0.9828 0.7614 0.5297 0.4533
## Proportion of Variance 0.4754 0.2853 0.1381 0.1012
## Cumulative Proportion  0.4754 0.7607 0.8988 1.0000
## 
## $item
## Importance of components:
##                          [,1]    [,2]
## Standard deviation     0.2824 0.08383
## Proportion of Variance 0.9190 0.08101
## Cumulative Proportion  0.9190 1.00000
# all components seem okay
### ANOVA model comparison - reduced 5th zcp - reduced 4th zcp model #############################
anova(zcp5_F02, zcp4_F02)
# not significant and AIC smaller for zcp5_F02 - preferred
# I go with zcp5_F02 as final zcp model and check correlations now
## REDUCED FINAL model with correlations =================================
summary(corr5_F02 <- lmer(f0_range2 ~ 
                            1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                            grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                            cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                            lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                            grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                            (1+grpd.ungr+rd.rpt+
                               grpd.ungr_cor.inc|Subj)+
                            (1+lhd.rhd|item),
                          REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## f0_range2 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + grpd.ungr + rd.rpt + grpd.ungr_cor.inc |  
##     Subj) + (1 + lhd.rhd | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   5934.6   6072.9  -2940.3   5880.6     1214 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0620 -0.4399 -0.0492  0.4312  6.7337 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr             
##  Subj     (Intercept)       3.42803  1.8515                    
##           grpd.ungr         1.11676  1.0568    0.57            
##           rd.rpt            1.66202  1.2892    0.05 -0.09      
##           grpd.ungr_cor.inc 5.26044  2.2936    0.57  0.14 -0.04
##  item     (Intercept)       0.04727  0.2174                    
##           lhd.rhd           0.37060  0.6088   -1.00            
##  Residual                   5.68054  2.3834                    
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                5.26272    0.34923  33.40963  15.069  < 2e-16 ***
## grpd.ungr                  2.05936    0.27595  29.61640   7.463 2.78e-08 ***
## lhd.rhd                    0.92840    0.64789  38.75062   1.433 0.159884    
## cor.inc                    0.39009    0.20503 793.88170   1.903 0.057458 .  
## rd.rpt                     2.37051    0.32570  37.53587   7.278 1.11e-08 ***
## grpd.ungr_lhd.rhd          0.04822    0.57323  32.83986   0.084 0.933469    
## grpd.ungr_cor.inc          3.23058    0.56081  29.55348   5.761 2.88e-06 ***
## cor.inc_lhd.rhd           -0.25374    0.40878 828.65554  -0.621 0.534951    
## grpd.ungr_rd.rpt           0.79556    0.36019  46.68496   2.209 0.032139 *  
## lhd.rhd_rd.rpt             2.82898    0.67214  39.57005   4.209 0.000143 ***
## cor.inc_rd.rpt             0.29084    0.38932 987.57499   0.747 0.455216    
## grpd.ungr_lhd.rhd_rd.rpt   1.41032    0.79387  31.45708   1.777 0.085324 .  
## cor.inc_lhd.rhd_rd.rpt     1.24264    0.77915 989.06389   1.595 0.111060    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
### ANOVA model comparison - reduced final model with correlation - initial model #############################
anova(corr5_F02, m0_F02)
# not significant and AIC smaller for corr5_F02 - preferred
# initial model does not accumulate more information
# FINAL model ---------------------------------
## model with factors and * for fixed effects =================================
# use REML=TRUE
summary(F02 <- lmer(f0_range2 ~ 
                      condition*accuracy+condition*group*exp+group*exp*accuracy+
                      (1+grpd.ungr+rd.rpt+
                         grpd.ungr_cor.inc|Subj)+
                      (1+lhd.rhd|item),
                    REML = TRUE,pDat_mm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: f0_range2 ~ condition * accuracy + condition * group * exp +  
##     group * exp * accuracy + (1 + grpd.ungr + rd.rpt + grpd.ungr_cor.inc |  
##     Subj) + (1 + lhd.rhd | item)
##    Data: pDat_mm
## 
## REML criterion at convergence: 5878.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0371 -0.4386 -0.0492  0.4286  6.7220 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr             
##  Subj     (Intercept)       3.65139  1.9109                    
##           grpd.ungr         1.27163  1.1277    0.55            
##           rd.rpt            1.83021  1.3529    0.05 -0.08      
##           grpd.ungr_cor.inc 5.58866  2.3640    0.56  0.11 -0.05
##  item     (Intercept)       0.05585  0.2363                    
##           lhd.rhd           0.43847  0.6622   -1.00            
##  Residual                   5.69990  2.3874                    
## Number of obs: 1241, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                         Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)              5.26412    0.36027  31.76563  14.612 1.19e-15 ***
## condition1               2.06366    0.28833  28.28438   7.157 8.17e-08 ***
## accuracy1                0.38055    0.20676 789.12333   1.841 0.066063 .  
## group1                   0.93450    0.67099  36.46869   1.393 0.172136    
## exp1                     2.37714    0.33821  34.85543   7.029 3.59e-08 ***
## condition1:accuracy1     3.20448    0.57272  28.60125   5.595 5.09e-06 ***
## condition1:group1        0.05207    0.60072  31.61980   0.087 0.931475    
## condition1:exp1          0.80868    0.36940  38.61847   2.189 0.034695 *  
## group1:exp1              2.84308    0.69995  37.01803   4.062 0.000243 ***
## accuracy1:group1        -0.25876    0.41212 821.92631  -0.628 0.530254    
## accuracy1:exp1           0.28597    0.39196 974.52750   0.730 0.465816    
## condition1:group1:exp1   1.42160    0.82438  26.74382   1.724 0.096167 .  
## accuracy1:group1:exp1    1.23295    0.78442 976.38331   1.572 0.116320    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# no warnings
save(F02, file = "F02.Rdata")
#*Nice table output*---------------------------------
# sjPlot::tab_model((F02), show.se = TRUE, show.stat = TRUE, collapse.ci=TRUE, file="F02.doc")
# clean workspace ---------------------------------
# ls()
rm(corr5_F02, m0_F02, zcp_F02, zcp1_F02, zcp2_F02, zcp3_F02, zcp4_F02, zcp5_F02)

4.4 Model reduction - pause (P)

# all "REML = FALSE" to use ANOVA model comparison later on
# INITIAL model with indicator variables ---------------------------------
summary(m0_P <- lmer(pause3 ~ 
                       1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                       grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                       cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                       lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                       grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                       (1+cor.inc+grpd.ungr+rd.rpt+
                          grpd.ungr_cor.inc+
                          cor.inc_rd.rpt+
                          grpd.ungr_rd.rpt|Subj)+
                       (1+lhd.rhd+cor.inc+
                          cor.inc_lhd.rhd|item),
                     REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## pause3 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + cor.inc_rd.rpt + grpd.ungr_rd.rpt | Subj) +  
##     (1 + lhd.rhd + cor.inc + cor.inc_lhd.rhd | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8029.1   8297.5  -3962.6   7925.1     1237 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5594 -0.3629 -0.0635  0.3633  6.6807 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)        7.2997  2.7018                                
##           cor.inc            4.3047  2.0748   -0.19                        
##           grpd.ungr          8.5614  2.9260    0.22  0.52                  
##           rd.rpt            12.8623  3.5864    0.27 -0.16  0.36            
##           grpd.ungr_cor.inc 50.7980  7.1273    0.78  0.07  0.09 -0.30      
##           cor.inc_rd.rpt    10.3345  3.2147    0.36  0.67  0.40 -0.46  0.76
##           grpd.ungr_rd.rpt  28.3623  5.3256    0.30  0.10  0.74  0.69  0.05
##  item     (Intercept)        0.2771  0.5264                                
##           lhd.rhd            0.1430  0.3782    0.08                        
##           cor.inc            0.6806  0.8250   -0.97  0.15                  
##           cor.inc_lhd.rhd    3.6581  1.9126    0.71 -0.64 -0.86            
##  Residual                   22.7076  4.7652                                
##       
##       
##       
##       
##       
##       
##       
##   0.07
##       
##       
##       
##       
##       
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                8.9587     0.5407 31.8062  16.568  < 2e-16 ***
## grpd.ungr                  8.2968     0.6928 31.1717  11.975 3.39e-13 ***
## lhd.rhd                    2.0314     0.8598 31.2238   2.363  0.02455 *  
## cor.inc                    2.0552     0.6091 30.2512   3.374  0.00204 ** 
## rd.rpt                     4.5842     0.8243 27.2811   5.561 6.54e-06 ***
## grpd.ungr_lhd.rhd          1.0749     1.3168 31.5545   0.816  0.42046    
## grpd.ungr_cor.inc         17.7254     1.5717 31.6610  11.278 1.27e-12 ***
## cor.inc_lhd.rhd           -0.2628     1.2298 28.8501  -0.214  0.83226    
## grpd.ungr_rd.rpt           7.3703     1.2044 27.8586   6.120 1.36e-06 ***
## lhd.rhd_rd.rpt            -0.7190     1.5242 28.8788  -0.472  0.64067    
## cor.inc_rd.rpt             0.5114     1.1437 33.1015   0.447  0.65769    
## grpd.ungr_lhd.rhd_rd.rpt  -0.6326     2.3589 27.7217  -0.268  0.79054    
## cor.inc_lhd.rhd_rd.rpt     2.0219     2.1337 36.4327   0.948  0.34956    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
# correlations seem okay for Subj and for items
summary(rePCA(m0_P))
## $Subj
## Importance of components:
##                          [,1]   [,2]   [,3]    [,4]    [,5]      [,6]      [,7]
## Standard deviation     1.6691 1.3482 0.7525 0.38927 0.27237 0.0005812 0.0001294
## Proportion of Variance 0.5163 0.3369 0.1050 0.02808 0.01375 0.0000000 0.0000000
## Cumulative Proportion  0.5163 0.8532 0.9582 0.98625 1.00000 1.0000000 1.0000000
## 
## $item
## Importance of components:
##                          [,1]    [,2]      [,3] [,4]
## Standard deviation     0.4398 0.12718 8.639e-05    0
## Proportion of Variance 0.9228 0.07718 0.000e+00    0
## Cumulative Proportion  0.9228 1.00000 1.000e+00    1
# 5 of 7 components for Subj | 2 of 4 components for item contribute (variance)
# try model reduction of randEff ---------------------------------
## start with ZERO CORRELATION model =================================
summary(zcp_P <- lmer(pause3 ~ 
                        1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                        grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                        cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                        lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                        grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                        (1+cor.inc+grpd.ungr+rd.rpt+
                           grpd.ungr_cor.inc+
                           cor.inc_rd.rpt+
                           grpd.ungr_rd.rpt||Subj)+
                        (1+lhd.rhd+cor.inc+
                           cor.inc_lhd.rhd||item),
                      REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## pause3 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + cor.inc_rd.rpt + grpd.ungr_rd.rpt ||  
##     Subj) + (1 + lhd.rhd + cor.inc + cor.inc_lhd.rhd || item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8063.9   8192.9  -4006.9   8013.9     1264 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4266 -0.3456 -0.0304  0.3188  6.8826 
## 
## Random effects:
##  Groups   Name              Variance  Std.Dev. 
##  Subj     (Intercept)       7.849e+00 2.8015381
##  Subj.1   cor.inc           4.978e+00 2.2311549
##  Subj.2   grpd.ungr         9.762e+00 3.1244283
##  Subj.3   rd.rpt            8.867e+00 2.9776854
##  Subj.4   grpd.ungr_cor.inc 4.453e+01 6.6730943
##  Subj.5   cor.inc_rd.rpt    7.430e-08 0.0002726
##  Subj.6   grpd.ungr_rd.rpt  2.165e+01 4.6529683
##  item     (Intercept)       0.000e+00 0.0000000
##  item.1   lhd.rhd           0.000e+00 0.0000000
##  item.2   cor.inc           0.000e+00 0.0000000
##  item.3   cor.inc_lhd.rhd   3.170e+00 1.7805591
##  Residual                   2.358e+01 4.8564295
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                9.0685     0.5564  31.7393  16.298  < 2e-16 ***
## grpd.ungr                  8.5559     0.7138  34.9168  11.986 6.27e-14 ***
## lhd.rhd                    2.5681     1.1139  31.8543   2.305 0.027808 *  
## cor.inc                    2.4926     0.6099  28.3101   4.087 0.000327 ***
## rd.rpt                     4.3801     0.7172  32.7562   6.107 7.25e-07 ***
## grpd.ungr_lhd.rhd          1.7753     1.4111  34.0911   1.258 0.216909    
## grpd.ungr_cor.inc         17.7410     1.4766  33.9227  12.015 8.99e-14 ***
## cor.inc_lhd.rhd            0.3076     1.2691  29.3077   0.242 0.810167    
## grpd.ungr_rd.rpt           6.9718     1.1236  28.4198   6.205 9.91e-07 ***
## lhd.rhd_rd.rpt            -1.0547     1.4312  32.4925  -0.737 0.466441    
## cor.inc_rd.rpt             0.1846     0.8994 624.1574   0.205 0.837460    
## grpd.ungr_lhd.rhd_rd.rpt  -0.6249     2.2559  29.1007  -0.277 0.783745    
## cor.inc_lhd.rhd_rd.rpt     1.3933     1.9449  88.5379   0.716 0.475638    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
    'Random effects:
     Groups   Name            Variance Std.Dev.
     Subj     (Intercept)         7.848   2.801   
     Subj.1   cor.inc             4.978   2.231   
     Subj.2   grpd.ungr           9.762   3.124   
     Subj.3   rd.rpt              8.867   2.978   
     Subj.4   grpd.ungr_cor.inc   44.529   6.673   
     Subj.5   cor.inc_rd.rpt      0.000   0.000 # no variance
     Subj.6   grpd.ungr_rd.rpt    21.651   4.653   
     item     (Intercept)         0.000   0.000   
     item.1   lhd.rhd             0.000   0.000 # no variance 
     item.2   cor.inc             0.000   0.000 # no variance
     item.3   cor.inc_lhd.rhd     3.171   1.781   
     Residual                 23.585   4.856   
    Number of obs: 1289, groups:  Subj, 32; item, 24'
## [1] "Random effects:\n     Groups   Name            Variance Std.Dev.\n     Subj     (Intercept)         7.848   2.801   \n     Subj.1   cor.inc             4.978   2.231   \n     Subj.2   grpd.ungr           9.762   3.124   \n     Subj.3   rd.rpt              8.867   2.978   \n     Subj.4   grpd.ungr_cor.inc   44.529   6.673   \n     Subj.5   cor.inc_rd.rpt      0.000   0.000 # no variance\n     Subj.6   grpd.ungr_rd.rpt    21.651   4.653   \n     item     (Intercept)         0.000   0.000   \n     item.1   lhd.rhd             0.000   0.000 # no variance \n     item.2   cor.inc             0.000   0.000 # no variance\n     item.3   cor.inc_lhd.rhd     3.171   1.781   \n     Residual                 23.585   4.856   \n    Number of obs: 1289, groups:  Subj, 32; item, 24"
## ANOVA model comparison - initial model - zcp model =================================
anova(m0_P, zcp_P)
# m0_P has signif. better fit - removing all correlations too harsh
## back to initial model =================================
# try exclusion of randEff slopes with regards to zcp
### FIRST #############################
# try exclusion of "cor.inc_rd.rpt" as slope per Subj
summary(m1_P <- lmer(pause3 ~ 
                       1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                       grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                       cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                       lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                       grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                       (1+cor.inc+grpd.ungr+rd.rpt+
                          grpd.ungr_cor.inc+
                          grpd.ungr_rd.rpt|Subj)+
                       (1+lhd.rhd+cor.inc+
                          cor.inc_lhd.rhd|item),
                     REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## pause3 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + grpd.ungr_rd.rpt | Subj) + (1 + lhd.rhd +  
##     cor.inc + cor.inc_lhd.rhd | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8024.4   8256.7  -3967.2   7934.4     1244 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5182 -0.3580 -0.0663  0.3493  6.6451 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)        7.80079 2.7930                                
##           cor.inc            5.09424 2.2570   -0.27                        
##           grpd.ungr          9.83907 3.1367    0.25  0.61                  
##           rd.rpt            11.06323 3.3261    0.39  0.01  0.44            
##           grpd.ungr_cor.inc 47.12526 6.8648    0.82 -0.11  0.15 -0.12      
##           grpd.ungr_rd.rpt  25.51835 5.0516    0.33 -0.01  0.65  0.80  0.02
##  item     (Intercept)        0.19577 0.4425                                
##           lhd.rhd            0.02583 0.1607   -0.41                        
##           cor.inc            0.58776 0.7667   -0.99  0.52                  
##           cor.inc_lhd.rhd    3.08416 1.7562    0.96 -0.66 -0.99            
##  Residual                   23.02807 4.7988                                
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)               8.949921   0.548281 31.719094  16.324  < 2e-16 ***
## grpd.ungr                 8.318558   0.713032 33.182192  11.666 2.76e-13 ***
## lhd.rhd                   1.949930   0.837642 33.498895   2.328   0.0261 *  
## cor.inc                   2.002837   0.626261 27.133272   3.198   0.0035 ** 
## rd.rpt                    4.823518   0.764060 33.070591   6.313 3.82e-07 ***
## grpd.ungr_lhd.rhd         1.375990   1.380263 32.592761   0.997   0.3262    
## grpd.ungr_cor.inc        17.611972   1.527892 33.443383  11.527 3.42e-13 ***
## cor.inc_lhd.rhd           0.157974   1.266159 25.962715   0.125   0.9017    
## grpd.ungr_rd.rpt          7.327576   1.169706 27.757274   6.264 9.36e-07 ***
## lhd.rhd_rd.rpt           -0.785435   1.459946 33.487118  -0.538   0.5941    
## cor.inc_rd.rpt           -0.007487   0.958254 89.913359  -0.008   0.9938    
## grpd.ungr_lhd.rhd_rd.rpt -0.486638   2.292581 28.667437  -0.212   0.8334    
## cor.inc_lhd.rhd_rd.rpt    2.563852   1.946562 68.057792   1.317   0.1922    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
summary(rePCA(m1_P))
## $Subj
## Importance of components:
##                          [,1]   [,2]    [,3]    [,4]    [,5] [,6]
## Standard deviation     1.5234 1.2950 0.65188 0.39631 0.20588    0
## Proportion of Variance 0.5021 0.3628 0.09194 0.03398 0.00917    0
## Cumulative Proportion  0.5021 0.8649 0.95685 0.99083 1.00000    1
## 
## $item
## Importance of components:
##                          [,1]    [,2]      [,3]      [,4]
## Standard deviation     0.4090 0.04284 8.107e-05 5.619e-06
## Proportion of Variance 0.9891 0.01085 0.000e+00 0.000e+00
## Cumulative Proportion  0.9891 1.00000 1.000e+00 1.000e+00
# 1 of 6 components for Subj | 2 of 4 components for item have zero variance
### ANOVA model comparison - initial model - reduced model #############################
anova(m0_P, m1_P)
# not significant and AIC smaller for m1_P - preferred
### SECOND #############################
# try exclusion also of "lhd.rhd" as slope per item
summary(m2_P <- lmer(pause3 ~ 
                       1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                       grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                       cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                       lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                       grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                       (1+cor.inc+grpd.ungr+rd.rpt+
                          grpd.ungr_cor.inc+
                          grpd.ungr_rd.rpt|Subj)+
                       (1+cor.inc+
                          cor.inc_lhd.rhd|item),
                     REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## pause3 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + grpd.ungr_rd.rpt | Subj) + (1 + cor.inc +  
##     cor.inc_lhd.rhd | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8016.4   8228.1  -3967.2   7934.4     1248 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5196 -0.3589 -0.0665  0.3494  6.6648 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)        7.8039  2.7935                                
##           cor.inc            5.0745  2.2527   -0.27                        
##           grpd.ungr          9.8340  3.1359    0.25  0.61                  
##           rd.rpt            11.0756  3.3280    0.39  0.01  0.44            
##           grpd.ungr_cor.inc 47.0909  6.8623    0.82 -0.11  0.15 -0.12      
##           grpd.ungr_rd.rpt  25.5208  5.0518    0.33  0.00  0.65  0.80  0.02
##  item     (Intercept)        0.2083  0.4565                                
##           cor.inc            0.6231  0.7893   -1.00                        
##           cor.inc_lhd.rhd    2.6435  1.6259    1.00 -1.00                  
##  Residual                   23.0385  4.7998                                
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                8.948821   0.548837  31.952675  16.305  < 2e-16 ***
## grpd.ungr                  8.310981   0.714354  35.258137  11.634 1.26e-13 ***
## lhd.rhd                    1.942741   0.834519  34.799515   2.328  0.02585 *  
## cor.inc                    2.001122   0.626867  27.475175   3.192  0.00352 ** 
## rd.rpt                     4.820944   0.765598  34.588301   6.297 3.30e-07 ***
## grpd.ungr_lhd.rhd          1.352954   1.378222  33.695812   0.982  0.33326    
## grpd.ungr_cor.inc         17.619918   1.528850  33.496345  11.525 3.36e-13 ***
## cor.inc_lhd.rhd            0.149601   1.257615  27.813327   0.119  0.90617    
## grpd.ungr_rd.rpt           7.344029   1.170218  28.759715   6.276 7.77e-07 ***
## lhd.rhd_rd.rpt            -0.775643   1.458619  34.119850  -0.532  0.59833    
## cor.inc_rd.rpt             0.005364   0.961115  93.109017   0.006  0.99556    
## grpd.ungr_lhd.rhd_rd.rpt  -0.417663   2.294091  29.058527  -0.182  0.85680    
## cor.inc_lhd.rhd_rd.rpt     2.573566   1.926683 123.584475   1.336  0.18409    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
# Model failed to converge with 1 negative eigenvalue: -1.7e+01
### THIRD #############################
# back to FIRST REDUCED model
# try exclusion also of "cor.inc" as slope per item
summary(m3_P <- lmer(pause3 ~ 
                       1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                       grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                       cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                       lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                       grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                       (1+cor.inc+grpd.ungr+rd.rpt+
                          grpd.ungr_cor.inc+
                          grpd.ungr_rd.rpt|Subj)+
                       (1+lhd.rhd+
                          cor.inc_lhd.rhd|item),
                     REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## pause3 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + grpd.ungr_rd.rpt | Subj) + (1 + lhd.rhd +  
##     cor.inc_lhd.rhd | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8019.0   8230.6  -3968.5   7937.0     1248 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5110 -0.3505 -0.0638  0.3542  6.8294 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)        7.76116 2.7859                                
##           cor.inc            5.37409 2.3182   -0.25                        
##           grpd.ungr          9.96203 3.1563    0.26  0.59                  
##           rd.rpt            10.84624 3.2934    0.40  0.01  0.45            
##           grpd.ungr_cor.inc 47.23947 6.8731    0.82 -0.09  0.16 -0.11      
##           grpd.ungr_rd.rpt  25.32193 5.0321    0.35 -0.05  0.65  0.80  0.02
##  item     (Intercept)        0.02525 0.1589                                
##           lhd.rhd            0.10694 0.3270   -0.84                        
##           cor.inc_lhd.rhd    4.65239 2.1569    0.94 -0.97                  
##  Residual                   23.17005 4.8135                                
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                8.9642     0.5409  30.6633  16.574  < 2e-16 ***
## grpd.ungr                  8.4258     0.6961  33.1977  12.104 1.01e-13 ***
## lhd.rhd                    1.9457     0.8511  34.9442   2.286  0.02843 *  
## cor.inc                    2.0548     0.6152  25.0731   3.340  0.00262 ** 
## rd.rpt                     4.7831     0.7410  33.6425   6.455 2.33e-07 ***
## grpd.ungr_lhd.rhd          1.4307     1.3963  33.2587   1.025  0.31291    
## grpd.ungr_cor.inc         17.5383     1.4858  31.2194  11.804 4.80e-13 ***
## cor.inc_lhd.rhd            0.2253     1.3092  26.0353   0.172  0.86470    
## grpd.ungr_rd.rpt           7.0596     1.1576  27.6135   6.099 1.49e-06 ***
## lhd.rhd_rd.rpt            -0.7651     1.4584  34.4065  -0.525  0.60321    
## cor.inc_rd.rpt            -0.1191     0.9099 345.4753  -0.131  0.89598    
## grpd.ungr_lhd.rhd_rd.rpt  -0.6383     2.2883  28.8280  -0.279  0.78227    
## cor.inc_lhd.rhd_rd.rpt     2.5114     2.0241  53.1778   1.241  0.22014    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
# Model failed to converge with max|grad| = 0.0207696 (tol = 0.002, component 1)
### FOURTH #############################
# back to FIRST REDUCED model
# try exclusion also of "cor.inc_lhd.rhd" as slope per item
summary(m4_P <- lmer(pause3 ~ 
                       1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                       grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                       cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                       lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                       grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                       (1+cor.inc+grpd.ungr+rd.rpt+
                          grpd.ungr_cor.inc+
                          grpd.ungr_rd.rpt|Subj)+
                       (1+lhd.rhd+cor.inc|item),
                     REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## pause3 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + grpd.ungr_rd.rpt | Subj) + (1 + lhd.rhd +  
##     cor.inc | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8019.5   8231.1  -3968.7   7937.5     1248 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4941 -0.3618 -0.0633  0.3446  6.8902 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)        7.6740  2.7702                                
##           cor.inc            5.0733  2.2524   -0.27                        
##           grpd.ungr          9.8985  3.1462    0.24  0.62                  
##           rd.rpt            11.2686  3.3569    0.39  0.01  0.43            
##           grpd.ungr_cor.inc 46.9611  6.8528    0.82 -0.11  0.15 -0.12      
##           grpd.ungr_rd.rpt  25.2664  5.0266    0.34  0.01  0.65  0.80  0.03
##  item     (Intercept)        0.3274  0.5722                                
##           lhd.rhd            0.2040  0.4517    1.00                        
##           cor.inc            0.9180  0.9581   -1.00 -1.00                  
##  Residual                   23.1687  4.8134                                
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                8.93784    0.55072  32.71151  16.229  < 2e-16 ***
## grpd.ungr                  8.28697    0.73114  34.69898  11.334 3.26e-13 ***
## lhd.rhd                    1.87398    0.82500  33.88069   2.271  0.02959 *  
## cor.inc                    2.03715    0.63883  27.47368   3.189  0.00355 ** 
## rd.rpt                     4.77708    0.78375  34.26532   6.095 6.29e-07 ***
## grpd.ungr_lhd.rhd          1.21257    1.39368  34.02248   0.870  0.39037    
## grpd.ungr_cor.inc         17.68798    1.53287  32.27654  11.539 5.43e-13 ***
## cor.inc_lhd.rhd            0.12183    1.21667  26.27744   0.100  0.92100    
## grpd.ungr_rd.rpt           7.32037    1.17762  28.78530   6.216 9.10e-07 ***
## lhd.rhd_rd.rpt            -0.68361    1.47888  34.35887  -0.462  0.64682    
## cor.inc_rd.rpt             0.08517    0.99386  58.35029   0.086  0.93200    
## grpd.ungr_lhd.rhd_rd.rpt  -0.27075    2.29671  29.12705  -0.118  0.90697    
## cor.inc_lhd.rhd_rd.rpt     2.66012    1.81970 355.05148   1.462  0.14467    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
summary(rePCA(m4_P))
## $Subj
## Importance of components:
##                          [,1]   [,2]    [,3]    [,4]    [,5]      [,6]
## Standard deviation     1.5147 1.2897 0.65459 0.38981 0.20764 9.061e-06
## Proportion of Variance 0.5008 0.3631 0.09353 0.03317 0.00941 0.000e+00
## Cumulative Proportion  0.5008 0.8639 0.95742 0.99059 1.00000 1.000e+00
## 
## $item
## Importance of components:
##                          [,1]      [,2] [,3]
## Standard deviation     0.2501 1.335e-06    0
## Proportion of Variance 1.0000 0.000e+00    0
## Cumulative Proportion  1.0000 1.000e+00    1
# 1 of 6 components for Subj | 2 of 3 components for item have zero variance
### ANOVA model comparison - initial model - reduced model #############################
anova(m0_P, m4_P)
# not significant and AIC smaller for m4_P - preferred
### FIFTH #############################
# try exclusion also of "lhd.rhd" as slope per item
summary(m5_P <- lmer(pause3 ~ 
                       1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                       grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                       cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                       lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                       grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                       (1+cor.inc+grpd.ungr+rd.rpt+
                          grpd.ungr_cor.inc+
                          grpd.ungr_rd.rpt|Subj)+
                       (1+cor.inc|item),
                     REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## pause3 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + grpd.ungr_rd.rpt | Subj) + (1 + cor.inc |      item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8015.0   8211.1  -3969.5   7939.0     1251 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4806 -0.3661 -0.0631  0.3485  6.8999 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)        7.6499  2.7658                                
##           cor.inc            5.0889  2.2558   -0.27                        
##           grpd.ungr          9.8426  3.1373    0.23  0.63                  
##           rd.rpt            11.1125  3.3335    0.39  0.01  0.44            
##           grpd.ungr_cor.inc 47.3394  6.8804    0.82 -0.11  0.15 -0.12      
##           grpd.ungr_rd.rpt  24.9721  4.9972    0.34  0.00  0.65  0.80  0.03
##  item     (Intercept)        0.1878  0.4333                                
##           cor.inc            0.7248  0.8514   -1.00                        
##  Residual                   23.2779  4.8247                                
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                8.93014    0.54435  31.68665  16.405  < 2e-16 ***
## grpd.ungr                  8.29449    0.71318  33.07472  11.630 3.14e-13 ***
## lhd.rhd                    1.86581    0.81542  33.58536   2.288  0.02855 *  
## cor.inc                    2.07624    0.63212  26.26237   3.285  0.00289 ** 
## rd.rpt                     4.75769    0.76472  32.81165   6.221 5.16e-07 ***
## grpd.ungr_lhd.rhd          1.24228    1.37846  33.47333   0.901  0.37391    
## grpd.ungr_cor.inc         17.72467    1.52533  31.54350  11.620 6.22e-13 ***
## cor.inc_lhd.rhd            0.15660    1.21608  25.91478   0.129  0.89853    
## grpd.ungr_rd.rpt           7.17001    1.15699  27.94072   6.197 1.09e-06 ***
## lhd.rhd_rd.rpt            -0.71160    1.46007  33.94266  -0.487  0.62913    
## cor.inc_rd.rpt             0.07505    0.97564  57.36875   0.077  0.93895    
## grpd.ungr_lhd.rhd_rd.rpt  -0.64562    2.27537  28.97303  -0.284  0.77862    
## cor.inc_lhd.rhd_rd.rpt     2.67822    1.81616 349.27549   1.475  0.14120    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
summary(rePCA(m5_P))
## $Subj
## Importance of components:
##                          [,1]   [,2]    [,3]    [,4]    [,5]      [,6]
## Standard deviation     1.5171 1.2790 0.65337 0.38784 0.19747 3.341e-06
## Proportion of Variance 0.5054 0.3592 0.09374 0.03303 0.00856 0.000e+00
## Cumulative Proportion  0.5054 0.8647 0.95841 0.99144 1.00000 1.000e+00
## 
## $item
## Importance of components:
##                         [,1]      [,2]
## Standard deviation     0.198 7.709e-06
## Proportion of Variance 1.000 0.000e+00
## Cumulative Proportion  1.000 1.000e+00
# 1 of 6 components for Subj | 1 of 2 components for item have zero variance
### ANOVA model comparison - m4_P - m5_P #############################
anova(m4_P, m5_P)
# not significant and AIC smaller for m5_P - preferred
### SIXTH #############################
# try exclusion also of "cor.inc" as slope per item - intercept only
summary(m6_P <- lmer(pause3 ~ 
                       1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                       grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                       cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                       lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                       grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                       (1+cor.inc+grpd.ungr+rd.rpt+
                          grpd.ungr_cor.inc+
                          grpd.ungr_rd.rpt|Subj)+
                       (1|item),
                     REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## pause3 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + grpd.ungr_rd.rpt | Subj) + (1 | item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8011.9   8197.8  -3970.0   7939.9     1253 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4598 -0.3677 -0.0655  0.3424  7.0755 
## 
## Random effects:
##  Groups   Name              Variance  Std.Dev. Corr                         
##  Subj     (Intercept)       7.702e+00 2.77532                               
##           cor.inc           5.089e+00 2.25584  -0.27                        
##           grpd.ungr         9.979e+00 3.15892   0.25  0.61                  
##           rd.rpt            1.094e+01 3.30787   0.41  0.02  0.46            
##           grpd.ungr_cor.inc 4.674e+01 6.83637   0.82 -0.12  0.15 -0.11      
##           grpd.ungr_rd.rpt  2.493e+01 4.99277   0.35 -0.01  0.66  0.81  0.03
##  item     (Intercept)       1.211e-08 0.00011                               
##  Residual                   2.347e+01 4.84438                               
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                8.95019    0.53759  30.70565  16.649  < 2e-16 ***
## grpd.ungr                  8.38150    0.69131  33.58594  12.124 8.13e-14 ***
## lhd.rhd                    1.89765    0.81713  34.17574   2.322  0.02630 *  
## cor.inc                    2.07156    0.60665  25.81942   3.415  0.00212 ** 
## rd.rpt                     4.74578    0.73918  33.92124   6.420 2.49e-07 ***
## grpd.ungr_lhd.rhd          1.31088    1.38334  33.55847   0.948  0.35010    
## grpd.ungr_cor.inc         17.62183    1.47667  31.20921  11.933 3.64e-13 ***
## cor.inc_lhd.rhd            0.16200    1.21488  25.91339   0.133  0.89495    
## grpd.ungr_rd.rpt           6.95256    1.14014  28.10458   6.098 1.39e-06 ***
## lhd.rhd_rd.rpt            -0.65959    1.45254  34.25120  -0.454  0.65262    
## cor.inc_rd.rpt            -0.03247    0.90747 351.22544  -0.036  0.97148    
## grpd.ungr_lhd.rhd_rd.rpt  -0.53740    2.27302  29.06302  -0.236  0.81476    
## cor.inc_lhd.rhd_rd.rpt     2.56241    1.81035 353.10480   1.415  0.15783    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
# Model failed to converge with 1 negative eigenvalue: -6.1e+00
### SEVENTH #############################
# BACK to m4_P
# try exclusion of "cor.inc" as slope per item instead of "lhd.rhd"
summary(m7_P <- lmer(pause3 ~ 
                       1+grpd.ungr+lhd.rhd+cor.inc+rd.rpt+
                       grpd.ungr_lhd.rhd+grpd.ungr_cor.inc+
                       cor.inc_lhd.rhd+grpd.ungr_rd.rpt+
                       lhd.rhd_rd.rpt+cor.inc_rd.rpt+
                       grpd.ungr_lhd.rhd_rd.rpt+cor.inc_lhd.rhd_rd.rpt+
                       (1+cor.inc+grpd.ungr+rd.rpt+
                          grpd.ungr_cor.inc+
                          grpd.ungr_rd.rpt|Subj)+
                       (1+lhd.rhd|item),
                     REML = FALSE,pDat_mm))
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## pause3 ~ 1 + grpd.ungr + lhd.rhd + cor.inc + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_cor.inc + cor.inc_lhd.rhd + grpd.ungr_rd.rpt +  
##     lhd.rhd_rd.rpt + cor.inc_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     cor.inc_lhd.rhd_rd.rpt + (1 + cor.inc + grpd.ungr + rd.rpt +  
##     grpd.ungr_cor.inc + grpd.ungr_rd.rpt | Subj) + (1 + lhd.rhd |      item)
##    Data: pDat_mm
## 
##      AIC      BIC   logLik deviance df.resid 
##   8015.7   8211.9  -3969.9   7939.7     1251 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4693 -0.3496 -0.0644  0.3397  7.0799 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)        7.71883 2.7783                                
##           cor.inc            5.07372 2.2525   -0.27                        
##           grpd.ungr         10.02589 3.1664    0.25  0.61                  
##           rd.rpt            10.99376 3.3157    0.40  0.01  0.45            
##           grpd.ungr_cor.inc 46.57779 6.8248    0.82 -0.12  0.15 -0.11      
##           grpd.ungr_rd.rpt  25.07541 5.0075    0.35 -0.02  0.66  0.81  0.03
##  item     (Intercept)        0.04957 0.2226                                
##           lhd.rhd            0.08933 0.2989   1.00                         
##  Residual                   23.40518 4.8379                                
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                8.95525    0.54003  30.91276  16.583  < 2e-16 ***
## grpd.ungr                  8.38203    0.69838  33.25103  12.002 1.24e-13 ***
## lhd.rhd                    1.90989    0.82162  34.10226   2.325  0.02619 *  
## cor.inc                    2.04901    0.60650  25.92744   3.378  0.00231 ** 
## rd.rpt                     4.75765    0.74604  33.68425   6.377 2.92e-07 ***
## grpd.ungr_lhd.rhd          1.30691    1.39098  33.54422   0.940  0.35416    
## grpd.ungr_cor.inc         17.59049    1.47540  31.22239  11.923 3.71e-13 ***
## cor.inc_lhd.rhd            0.13121    1.21458  26.04676   0.108  0.91480    
## grpd.ungr_rd.rpt           6.97248    1.15675  27.47555   6.028 1.83e-06 ***
## lhd.rhd_rd.rpt            -0.64188    1.45993  34.24185  -0.440  0.66294    
## cor.inc_rd.rpt            -0.04458    0.90792 353.73071  -0.049  0.96086    
## grpd.ungr_lhd.rhd_rd.rpt  -0.51727    2.29016  29.03852  -0.226  0.82289    
## cor.inc_lhd.rhd_rd.rpt     2.52673    1.81156 355.09139   1.395  0.16395    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
summary(rePCA(m7_P))
## $Subj
## Importance of components:
##                          [,1]   [,2]    [,3]    [,4]    [,5] [,6]
## Standard deviation     1.5044 1.2790 0.64670 0.38973 0.19191    0
## Proportion of Variance 0.5023 0.3630 0.09281 0.03371 0.00817    0
## Cumulative Proportion  0.5023 0.8653 0.95812 0.99183 1.00000    1
## 
## $item
## Importance of components:
##                           [,1]      [,2]
## Standard deviation     0.07704 1.904e-05
## Proportion of Variance 1.00000 0.000e+00
## Cumulative Proportion  1.00000 1.000e+00
# 1 of 6 components for Subj | 1 of 2 components for item have zero variance
# same for m5_P
### ANOVA model comparison - initial model - m7_P #############################
anova(m0_P, m7_P)
# not significant and AIC smaller for m7_P - preferred
### ANOVA model comparison - initial model - m5_P #############################
anova(m0_P, m5_P)
# not significant and AIC smaller for m7_P - preferred
# both ANOVAs obtain almost identical measures
# stay with m7_P
# FINAL model ---------------------------------
## model with factors and * for fixed effects =================================
# use REML=TRUE
summary(P <- lmer(pause3 ~ 
                    condition*accuracy+condition*group*exp+
                    group*exp*accuracy+
                    (1+cor.inc+grpd.ungr+rd.rpt+
                       grpd.ungr_cor.inc+
                       grpd.ungr_rd.rpt|Subj)+
                    (1+lhd.rhd|item),
                  REML = TRUE,pDat_mm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pause3 ~ condition * accuracy + condition * group * exp + group *  
##     exp * accuracy + (1 + cor.inc + grpd.ungr + rd.rpt + grpd.ungr_cor.inc +  
##     grpd.ungr_rd.rpt | Subj) + (1 + lhd.rhd | item)
##    Data: pDat_mm
## 
## REML criterion at convergence: 7918.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4809 -0.3452 -0.0567  0.3338  7.0880 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr                         
##  Subj     (Intercept)        8.25955 2.8739                                
##           cor.inc            6.02594 2.4548   -0.27                        
##           grpd.ungr         11.16365 3.3412    0.26  0.57                  
##           rd.rpt            12.10076 3.4786    0.40  0.01  0.44            
##           grpd.ungr_cor.inc 48.73116 6.9808    0.81 -0.10  0.14 -0.12      
##           grpd.ungr_rd.rpt  28.19503 5.3099    0.34 -0.06  0.63  0.80  0.01
##  item     (Intercept)        0.08123 0.2850                                
##           lhd.rhd            0.14436 0.3799   1.00                         
##  Residual                   23.41247 4.8386                                
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                        Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)              8.9745     0.5588  29.7961  16.060 3.25e-16 ***
## condition1               8.4361     0.7327  31.6435  11.514 7.54e-13 ***
## accuracy1                2.0183     0.6377  23.9710   3.165  0.00418 ** 
## group1                   1.9202     0.8571  31.7202   2.240  0.03220 *  
## exp1                     4.8107     0.7782  31.9286   6.182 6.48e-07 ***
## condition1:accuracy1    17.4744     1.5081  29.9484  11.587 1.36e-12 ***
## condition1:group1        1.3088     1.4557  31.5896   0.899  0.37541    
## condition1:exp1          6.9831     1.2164  26.3545   5.741 4.58e-06 ***
## group1:exp1             -0.5902     1.5188  31.9723  -0.389  0.70016    
## accuracy1:group1         0.1183     1.2768  24.2112   0.093  0.92696    
## accuracy1:exp1          -0.1541     0.9227 362.3933  -0.167  0.86743    
## condition1:group1:exp1  -0.4799     2.3977  27.5151  -0.200  0.84282    
## accuracy1:group1:exp1    2.4604     1.8416 362.2902   1.336  0.18238    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
save(P, file = "P.Rdata")
# *Nice table output* ---------------------------------
# sjPlot::tab_model((P), show.se = TRUE, show.stat = TRUE, collapse.ci=TRUE, file="P.doc")
# clean workspace ---------------------------------
# ls()
rm(m0_P, m1_P, m2_P, m3_P, m4_P, m5_P, m6_P, m7_P, zcp_P)

5 Final lengthening model (FL2) for name2

5.1 Visualizations and Post-Hoc-Tests

# interaction visualization ---------------------------------
# condition:accuracy =================================
### condition ~ accuracy #############################
load("FL2.Rdata")
emmeans::emmip(FL2, condition ~ accuracy, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - % final lengthening on name2 \n(relative to duration of name2)",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

### accuracy ~ condition #############################
emmeans::emmip(FL2, accuracy ~ condition, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - % final lengthening on name2 \n(relative to duration of name2)",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

# POST-HOC testing ---------------------------------
# what are the significant differences
emmeans::emmeans(FL2, pairwise ~ condition*accuracy) %>% 
    summary(by = NULL, adjust = "bonferroni")
## $emmeans
##  condition accuracy  emmean   SE   df lower.CL upper.CL
##  grouped   correct     48.1 1.21 35.7     44.9     51.3
##  ungrouped correct     38.7 1.48 49.6     34.8     42.5
##  grouped   incorrect   43.5 1.19 36.1     40.4     46.6
##  ungrouped incorrect   43.5 1.40 37.4     39.9     47.2
## 
## Results are averaged over the levels of: group, exp 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## 
## $contrasts
##  contrast                                estimate   SE   df t.ratio p.value
##  grouped correct - ungrouped correct       9.4548 1.69 39.3   5.607  <.0001
##  grouped correct - grouped incorrect       4.6467 1.02 29.4   4.547  0.0005
##  grouped correct - ungrouped incorrect     4.5924 1.53 42.1   2.993  0.0276
##  ungrouped correct - grouped incorrect    -4.8081 1.49 38.0  -3.227  0.0155
##  ungrouped correct - ungrouped incorrect  -4.8624 1.23 34.1  -3.948  0.0022
##  grouped incorrect - ungrouped incorrect  -0.0543 1.46 32.0  -0.037  1.0000
## 
## Results are averaged over the levels of: group, exp 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: bonferroni method for 6 tests
# interaction visualization ---------------------------------
# group:exp =================================
### group ~ exp #############################
emmeans::emmip(FL2, group ~ exp, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - % final lengthening on name2 \n(relative to duration of name2)",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

### exp ~ group #############################
emmeans::emmip(FL2, exp ~ group, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - % final lengthening on name2 \n(relative to duration of name2)",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

# POST-HOC testing ---------------------------------
# what are the significant differences
emmeans::emmeans(FL2, pairwise ~ group*exp) %>% 
    summary(by = NULL, adjust = "bonferroni")
## $emmeans
##  group exp           emmean   SE   df lower.CL upper.CL
##  LHDP  reading_aloud   45.3 1.49 37.7     41.3     49.2
##  RHDP  reading_aloud   43.4 1.67 44.2     39.1     47.7
##  LHDP  repetition      45.7 1.54 34.7     41.6     49.7
##  RHDP  repetition      39.6 1.65 42.3     35.2     43.9
## 
## Results are averaged over the levels of: condition, accuracy 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## 
## $contrasts
##  contrast                                estimate   SE   df t.ratio p.value
##  LHDP reading_aloud - RHDP reading_aloud    1.852 1.99 37.1   0.932  1.0000
##  LHDP reading_aloud - LHDP repetition      -0.405 1.51 38.3  -0.268  1.0000
##  LHDP reading_aloud - RHDP repetition       5.700 2.25 58.6   2.538  0.0829
##  RHDP reading_aloud - LHDP repetition      -2.257 2.29 59.4  -0.987  1.0000
##  RHDP reading_aloud - RHDP repetition       3.848 1.90 28.6   2.021  0.3163
##  LHDP repetition - RHDP repetition          6.104 2.03 34.6   3.013  0.0289
## 
## Results are averaged over the levels of: condition, accuracy 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: bonferroni method for 6 tests

6 F0_range model (F02) for name2

6.1 Visualizations and Post-Hoc-Tests

# interaction visualization ---------------------------------
load("F02.Rdata")
# condition:accuracy =================================
### condition ~ accuracy #############################
emmeans::emmip(F02, condition ~ accuracy, 
               CIs = TRUE, CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - f0 range on name2 in semitones",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

### accuracy ~ condition #############################
emmeans::emmip(F02, accuracy ~ condition, 
               CIs = TRUE, CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - f0 range on name2 in semitones",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

# POST-HOC testing ---------------------------------
# what are the significant differences
emmeans::emmeans(F02, pairwise ~ condition*accuracy) %>% 
    summary(by = NULL, adjust = "bonferroni")
## $emmeans
##  condition accuracy  emmean    SE   df lower.CL upper.CL
##  grouped   correct     7.29 0.500 33.1     5.97     8.61
##  ungrouped correct     3.62 0.289 32.9     2.86     4.38
##  grouped   incorrect   5.30 0.452 38.8     4.12     6.49
##  ungrouped incorrect   4.84 0.451 39.4     3.66     6.02
## 
## Results are averaged over the levels of: group, exp 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## 
## $contrasts
##  contrast                                estimate    SE   df t.ratio p.value
##  grouped correct - ungrouped correct        3.666 0.374 33.7   9.790  <.0001
##  grouped correct - grouped incorrect        1.983 0.380 71.7   5.224  <.0001
##  grouped correct - ungrouped incorrect      2.444 0.370 75.2   6.601  <.0001
##  ungrouped correct - grouped incorrect     -1.683 0.355 51.6  -4.743  0.0001
##  ungrouped correct - ungrouped incorrect   -1.222 0.343 58.6  -3.562  0.0044
##  grouped incorrect - ungrouped incorrect    0.461 0.449 28.6   1.028  1.0000
## 
## Results are averaged over the levels of: group, exp 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: bonferroni method for 6 tests
# interaction visualization ---------------------------------
# condition:exp =================================
### condition ~ exp #############################
emmeans::emmip(F02, condition ~ exp, 
               CIs = TRUE, CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - f0 range on name2 in semitones",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

### exp ~ condition #############################
emmeans::emmip(F02, exp ~ condition, 
               CIs = TRUE, CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - f0 range on name2 in semitones",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

# POST-HOC testing ---------------------------------
# what are the significant differences
emmeans::emmeans(F02, pairwise ~ condition*exp) %>% 
    summary(by = NULL, adjust = "bonferroni")
## $emmeans
##  condition exp           emmean    SE   df lower.CL upper.CL
##  grouped   reading_aloud   7.69 0.479 36.6     6.43     8.95
##  ungrouped reading_aloud   5.22 0.389 38.7     4.20     6.24
##  grouped   repetition      4.91 0.481 33.1     3.63     6.18
##  ungrouped repetition      3.25 0.390 34.2     2.22     4.27
## 
## Results are averaged over the levels of: accuracy, group 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## 
## $contrasts
##  contrast                                        estimate    SE   df t.ratio
##  grouped reading_aloud - ungrouped reading_aloud    2.468 0.341 36.6   7.227
##  grouped reading_aloud - grouped repetition         2.781 0.396 43.5   7.016
##  grouped reading_aloud - ungrouped repetition       4.441 0.451 33.4   9.854
##  ungrouped reading_aloud - grouped repetition       0.313 0.454 32.0   0.690
##  ungrouped reading_aloud - ungrouped repetition     1.973 0.388 42.4   5.082
##  grouped repetition - ungrouped repetition          1.659 0.353 40.4   4.696
##  p.value
##   <.0001
##   <.0001
##   <.0001
##   1.0000
##   <.0001
##   0.0002
## 
## Results are averaged over the levels of: accuracy, group 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: bonferroni method for 6 tests
# interaction visualization ---------------------------------
# group:exp =================================
### exp ~ group #############################
emmeans::emmip(F02, exp ~ group, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - f0 range on name2 in semitones",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

### group ~ exp #############################
emmeans::emmip(F02, group ~ exp, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - f0 range on name2 in semitones",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

# POST-HOC testing ---------------------------------
# what are the significant differences
emmeans::emmeans(F02, pairwise ~ group*exp) %>% 
    summary(by = NULL, adjust = "bonferroni")
## $emmeans
##  group exp           emmean    SE   df lower.CL upper.CL
##  LHDP  reading_aloud   7.63 0.597 38.7     6.07     9.20
##  RHDP  reading_aloud   5.27 0.535 39.2     3.87     6.67
##  LHDP  repetition      3.83 0.593 32.3     2.26     5.40
##  RHDP  repetition      4.32 0.519 34.8     2.95     5.69
## 
## Results are averaged over the levels of: condition, accuracy 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## 
## $contrasts
##  contrast                                estimate    SE   df t.ratio p.value
##  LHDP reading_aloud - RHDP reading_aloud    2.356 0.800 39.1   2.944  0.0326
##  LHDP reading_aloud - LHDP repetition       3.799 0.532 34.3   7.140  <.0001
##  LHDP reading_aloud - RHDP repetition       3.312 0.778 51.6   4.258  0.0005
##  RHDP reading_aloud - LHDP repetition       1.443 0.779 47.8   1.851  0.4218
##  RHDP reading_aloud - RHDP repetition       0.956 0.470 34.3   2.032  0.2998
##  LHDP repetition - RHDP repetition         -0.487 0.775 31.6  -0.629  1.0000
## 
## Results are averaged over the levels of: condition, accuracy 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: bonferroni method for 6 tests

7 Pause model (P) for pause after name2

7.1 Visualizations and Post-Hoc-Tests

# interaction visualization ---------------------------------
load("P.Rdata")
# condition:accuracy =================================
### condition ~ accuracy #############################
emmeans::emmip(P, condition ~ accuracy, 
               CIs = TRUE, CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - % pause after name2 \n(relative to utterance duration)",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

### accuracy ~ condition #############################
emmeans::emmip(P, accuracy ~ condition, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - % pause after name2 \n(relative to utterance duration)",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

# POST-HOC testing ---------------------------------
# what are the significant differences
emmeans::emmeans(P, pairwise ~ condition * accuracy) %>% 
    summary(by = NULL, adjust = "bonferroni")
## $emmeans
##  condition accuracy  emmean    SE   df lower.CL upper.CL
##  grouped   correct    18.57 1.036 29.4    15.81    21.33
##  ungrouped correct     1.40 0.358 27.0     0.44     2.35
##  grouped   incorrect   7.81 0.774 31.1     5.76     9.87
##  ungrouped incorrect   8.12 1.056 28.9     5.30    10.93
## 
## Results are averaged over the levels of: group, exp 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## 
## $contrasts
##  contrast                                estimate    SE   df t.ratio p.value
##  grouped correct - ungrouped correct       17.173 1.055 29.6  16.278  <.0001
##  grouped correct - grouped incorrect       10.755 1.053 28.8  10.217  <.0001
##  grouped correct - ungrouped incorrect     10.454 1.183 28.5   8.835  <.0001
##  ungrouped correct - grouped incorrect     -6.418 0.802 31.2  -7.999  <.0001
##  ungrouped correct - ungrouped incorrect   -6.719 0.999 28.5  -6.726  <.0001
##  grouped incorrect - ungrouped incorrect   -0.301 1.116 30.3  -0.270  1.0000
## 
## Results are averaged over the levels of: group, exp 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: bonferroni method for 6 tests
# interaction visualization ---------------------------------
# condition:exp =================================
### condition ~ exp #############################
emmeans::emmip(P, condition ~ exp, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - % pause after name2 \n(relative to utterance duration)",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

### exp ~ condition #############################
emmeans::emmip(P, exp ~ condition, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - % pause after name2 \n(relative to utterance duration)",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

# POST-HOC testing ---------------------------------
# what are the significant differences
emmeans::emmeans(P, pairwise ~ condition * exp) %>% 
    summary(by = NULL, adjust = "bonferroni")
## $emmeans
##  condition exp           emmean    SE   df lower.CL upper.CL
##  grouped   reading_aloud  17.34 1.165 31.7    14.26    20.43
##  ungrouped reading_aloud   5.42 0.729 32.3     3.49     7.34
##  grouped   repetition      9.04 0.744 27.3     7.05    11.03
##  ungrouped repetition      4.10 0.674 30.7     2.31     5.88
## 
## Results are averaged over the levels of: accuracy, group 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## 
## $contrasts
##  contrast                                        estimate    SE   df t.ratio
##  grouped reading_aloud - ungrouped reading_aloud    11.93 1.131 29.6  10.547
##  grouped reading_aloud - grouped repetition          8.30 1.258 29.4   6.601
##  grouped reading_aloud - ungrouped repetition       13.25 1.275 35.4  10.390
##  ungrouped reading_aloud - grouped repetition       -3.63 0.897 29.7  -4.041
##  ungrouped reading_aloud - ungrouped repetition      1.32 0.695 38.0   1.899
##  grouped repetition - ungrouped repetition           4.94 0.807 26.2   6.124
##  p.value
##   <.0001
##   <.0001
##   <.0001
##   0.0021
##   0.3911
##   <.0001
## 
## Results are averaged over the levels of: accuracy, group 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: bonferroni method for 6 tests

8 Model reduction - production accuracy

# INITIAL MAX model ---------------------------------
summary(rateM_i <- glmer(rate ~ condition*group*exp+
                           (1+condition*exp|Subj)+
                           (1+group|item),
                         family="binomial",
                         control=glmerControl(optCtrl=list(maxfun=1e5),
                                              optimizer="bobyqa"),
                         pDat))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: rate ~ condition * group * exp + (1 + condition * exp | Subj) +  
##     (1 + group | item)
##    Data: pDat
## Control: glmerControl(optCtrl = list(maxfun = 1e+05), optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1172.7   1281.1   -565.4   1130.7     1268 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5739 -0.1593  0.2867  0.4828  3.0638 
## 
## Random effects:
##  Groups Name            Variance Std.Dev. Corr           
##  Subj   (Intercept)     2.57116  1.6035                  
##         condition1      8.23428  2.8695   0.95           
##         exp1            4.13984  2.0347   0.34  0.38     
##         condition1:exp1 4.38362  2.0937   0.11  0.23 0.95
##  item   (Intercept)     0.01638  0.1280                  
##         group1          0.13368  0.3656   -1.00          
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              1.5818     0.3629   4.359 1.31e-05 ***
## condition1               0.1798     0.6630   0.271 0.786291    
## group1                   0.7818     0.6713   1.165 0.244157    
## exp1                     2.3320     0.5704   4.089 4.34e-05 ***
## condition1:group1        1.1029     1.2059   0.915 0.360408    
## condition1:exp1          3.5188     0.9340   3.767 0.000165 ***
## group1:exp1              1.4909     1.0070   1.481 0.138704    
## condition1:group1:exp1   2.2744     1.5297   1.487 0.137071    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtn1 group1 exp1  cndtn1:g1 cndtn1:x1 grp1:1
## condition1  0.902                                                
## group1      0.248  0.239                                         
## exp1        0.426  0.453  0.143                                  
## cndtn1:grp1 0.242  0.256  0.887  0.165                           
## condtn1:xp1 0.356  0.419  0.177  0.823 0.197                     
## group1:exp1 0.148  0.161  0.281  0.351 0.319     0.366           
## cndtn1:g1:1 0.188  0.200  0.167  0.395 0.225     0.476     0.771 
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# 'isSingular'
# all "REML = FALSE" to use ANOVA model comparison later on
# INITIAL model with indicator variables ---------------------------------
summary(m0_rateM <- glmer(rate ~ 1+
                            grpd.ungr+lhd.rhd+rd.rpt+
                            grpd.ungr_lhd.rhd+grpd.ungr_rd.rpt+lhd.rhd_rd.rpt+
                            grpd.ungr_lhd.rhd_rd.rpt+
                            (1+grpd.ungr+rd.rpt+grpd.ungr_rd.rpt|Subj)+
                            (1+lhd.rhd|item),
                          family="binomial",
                          control=glmerControl(optCtrl=list(maxfun=1e5),
                                               optimizer="bobyqa"),
                          pDat_mm))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: rate ~ 1 + grpd.ungr + lhd.rhd + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_rd.rpt + lhd.rhd_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     (1 + grpd.ungr + rd.rpt + grpd.ungr_rd.rpt | Subj) + (1 +  
##     lhd.rhd | item)
##    Data: pDat_mm
## Control: glmerControl(optCtrl = list(maxfun = 1e+05), optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1172.7   1281.1   -565.4   1130.7     1268 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5739 -0.1593  0.2867  0.4828  3.0638 
## 
## Random effects:
##  Groups Name             Variance Std.Dev. Corr           
##  Subj   (Intercept)      2.57116  1.6035                  
##         grpd.ungr        8.23428  2.8695   0.95           
##         rd.rpt           4.13984  2.0347   0.34  0.38     
##         grpd.ungr_rd.rpt 4.38362  2.0937   0.11  0.23 0.95
##  item   (Intercept)      0.01638  0.1280                  
##         lhd.rhd          0.13368  0.3656   -1.00          
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.5818     0.3629   4.359 1.31e-05 ***
## grpd.ungr                  0.1798     0.6630   0.271 0.786291    
## lhd.rhd                    0.7818     0.6713   1.165 0.244157    
## rd.rpt                     2.3320     0.5704   4.089 4.34e-05 ***
## grpd.ungr_lhd.rhd          1.1029     1.2059   0.915 0.360408    
## grpd.ungr_rd.rpt           3.5188     0.9340   3.767 0.000165 ***
## lhd.rhd_rd.rpt             1.4909     1.0070   1.481 0.138704    
## grpd.ungr_lhd.rhd_rd.rpt   2.2744     1.5297   1.487 0.137071    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grpd.n lhd.rh rd.rpt grpd.ngr_l. grpd.ngr_r. lhd._.
## grpd.ungr   0.902                                                     
## lhd.rhd     0.248  0.239                                              
## rd.rpt      0.426  0.453  0.143                                       
## grpd.ngr_l. 0.242  0.256  0.887  0.165                                
## grpd.ngr_r. 0.356  0.419  0.177  0.823  0.197                         
## lhd.rhd_rd. 0.148  0.161  0.281  0.351  0.319       0.366             
## grpd.ng_._. 0.188  0.200  0.167  0.395  0.225       0.476       0.771 
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# same model output
# two correlations close to 1 for Subj | -1 for items
summary(rePCA(m0_rateM))
## $Subj
## Importance of components:
##                          [,1]   [,2]    [,3] [,4]
## Standard deviation     3.5300 2.5583 0.56875    0
## Proportion of Variance 0.6447 0.3386 0.01674    0
## Cumulative Proportion  0.6447 0.9833 1.00000    1
## 
## $item
## Importance of components:
##                          [,1]      [,2]
## Standard deviation     0.3874 3.472e-07
## Proportion of Variance 1.0000 0.000e+00
## Cumulative Proportion  1.0000 1.000e+00
# 3 of 4 components for Subj | only intercept for item contribute (variance)
# try model reduction of randEff ---------------------------------
## start with ZERO CORRELATION model =================================
summary(zcp_rateM <- glmer(rate ~ 1+
                             grpd.ungr+lhd.rhd+rd.rpt+
                             grpd.ungr_lhd.rhd+grpd.ungr_rd.rpt+lhd.rhd_rd.rpt+
                             grpd.ungr_lhd.rhd_rd.rpt+
                             (1+grpd.ungr+rd.rpt+grpd.ungr_rd.rpt||Subj)+
                             (1+lhd.rhd||item),
                           family="binomial",
                           control=glmerControl(optCtrl=list(maxfun=1e5),
                                                optimizer="bobyqa"),
                           pDat_mm))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: rate ~ 1 + grpd.ungr + lhd.rhd + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_rd.rpt + lhd.rhd_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     (1 + grpd.ungr + rd.rpt + grpd.ungr_rd.rpt || Subj) + (1 +  
##     lhd.rhd || item)
##    Data: pDat_mm
## Control: glmerControl(optCtrl = list(maxfun = 1e+05), optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1202.9   1275.2   -587.5   1174.9     1275 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2327 -0.2725  0.2721  0.4679  2.3564 
## 
## Random effects:
##  Groups Name             Variance Std.Dev.
##  Subj   (Intercept)      1.57842  1.2564  
##  Subj.1 grpd.ungr        4.03538  2.0088  
##  Subj.2 rd.rpt           2.28673  1.5122  
##  Subj.3 grpd.ungr_rd.rpt 1.75237  1.3238  
##  item   (Intercept)      0.01059  0.1029  
##  item.1 lhd.rhd          0.18332  0.4282  
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.4294     0.2607   5.484 4.16e-08 ***
## grpd.ungr                 -0.5452     0.4223  -1.291    0.197    
## lhd.rhd                    0.5521     0.5202   1.061    0.289    
## rd.rpt                     1.6939     0.3808   4.448 8.66e-06 ***
## grpd.ungr_lhd.rhd          0.7744     0.8570   0.904    0.366    
## grpd.ungr_rd.rpt           2.1808     0.5363   4.066 4.78e-05 ***
## lhd.rhd_rd.rpt             1.1868     0.7656   1.550    0.121    
## grpd.ungr_lhd.rhd_rd.rpt   0.8611     1.0715   0.804    0.422    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grpd.n lhd.rh rd.rpt grpd.ngr_l. grpd.ngr_r. lhd._.
## grpd.ungr    0.045                                                    
## lhd.rhd      0.187  0.028                                             
## rd.rpt       0.014  0.049 -0.020                                      
## grpd.ngr_l.  0.036  0.184  0.041  0.049                               
## grpd.ngr_r.  0.098  0.032  0.051  0.175  0.031                        
## lhd.rhd_rd.  0.000  0.043 -0.033  0.289  0.047       0.118            
## grpd.ng_._.  0.056  0.019  0.047  0.100  0.018       0.350       0.124
    'Random effects:
     Groups Name           Variance Std.Dev.
     Subj   (Intercept)       1.57842  1.2564  
     Subj.1 grpd.ungr         4.03538  2.0088  
     Subj.2 rd.rpt            2.28673  1.5122  
     Subj.3 grpd.ungr_rd.rpt  1.75237  1.3238 # mean var. almost equal to sd 
     item   (Intercept)       0.01059  0.1029 # almost no variance 
     item.1 lhd.rhd           0.18332  0.4282 # sd larger than variance  
    Number of obs: 1289, groups:  Subj, 32; item, 24'
## [1] "Random effects:\n     Groups Name           Variance Std.Dev.\n     Subj   (Intercept)       1.57842  1.2564  \n     Subj.1 grpd.ungr         4.03538  2.0088  \n     Subj.2 rd.rpt            2.28673  1.5122  \n     Subj.3 grpd.ungr_rd.rpt  1.75237  1.3238 # mean var. almost equal to sd \n     item   (Intercept)       0.01059  0.1029 # almost no variance \n     item.1 lhd.rhd           0.18332  0.4282 # sd larger than variance  \n    Number of obs: 1289, groups:  Subj, 32; item, 24"
## ANOVA model comparison - initial model - zcp model =================================
anova(m0_rateM, zcp_rateM)
# m0_rateM has signif. better fit - removing all correlations too harsh
## back to initial model =================================
# try exclusion of randEff slopes with regards to zcp
### FIRST #############################
# try exclusion of "lhd.rhd" as slope per item
summary(m1_rateM <- glmer(rate ~ 1+
                            grpd.ungr+lhd.rhd+rd.rpt+
                            grpd.ungr_lhd.rhd+grpd.ungr_rd.rpt+lhd.rhd_rd.rpt+
                            grpd.ungr_lhd.rhd_rd.rpt+
                            (1+grpd.ungr+rd.rpt+grpd.ungr_rd.rpt|Subj)+
                            (1|item),
                          family="binomial",
                          control=glmerControl(optCtrl=list(maxfun=1e5),
                                               optimizer="bobyqa"),
                          pDat_mm))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: rate ~ 1 + grpd.ungr + lhd.rhd + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_rd.rpt + lhd.rhd_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     (1 + grpd.ungr + rd.rpt + grpd.ungr_rd.rpt | Subj) + (1 |      item)
##    Data: pDat_mm
## Control: glmerControl(optCtrl = list(maxfun = 1e+05), optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1170.2   1268.3   -566.1   1132.2     1270 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6493 -0.1576  0.2865  0.4938  3.1169 
## 
## Random effects:
##  Groups Name             Variance Std.Dev. Corr          
##  Subj   (Intercept)      2.53998  1.5937                 
##         grpd.ungr        8.12449  2.8503   0.95          
##         rd.rpt           4.06650  2.0166   0.34 0.38     
##         grpd.ungr_rd.rpt 4.31282  2.0767   0.12 0.24 0.95
##  item   (Intercept)      0.01979  0.1407                 
## Number of obs: 1289, groups:  Subj, 32; item, 24
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.5740     0.3613   4.356 1.32e-05 ***
## grpd.ungr                  0.1844     0.6599   0.279 0.779904    
## lhd.rhd                    0.7963     0.6634   1.200 0.229971    
## rd.rpt                     2.3249     0.5680   4.093 4.26e-05 ***
## grpd.ungr_lhd.rhd          1.0842     1.1897   0.911 0.362133    
## grpd.ungr_rd.rpt           3.5059     0.9319   3.762 0.000168 ***
## lhd.rhd_rd.rpt             1.5060     0.9897   1.522 0.128083    
## grpd.ungr_lhd.rhd_rd.rpt   2.2923     1.4935   1.535 0.124822    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grpd.n lhd.rh rd.rpt grpd.ngr_l. grpd.ngr_r. lhd._.
## grpd.ungr   0.901                                                     
## lhd.rhd     0.260  0.241                                              
## rd.rpt      0.430  0.456  0.147                                       
## grpd.ngr_l. 0.244  0.270  0.900  0.169                                
## grpd.ngr_r. 0.359  0.421  0.181  0.822  0.200                         
## lhd.rhd_rd. 0.153  0.166  0.291  0.372  0.329       0.373             
## grpd.ng_._. 0.195  0.206  0.176  0.405  0.235       0.510       0.796 
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
summary(rePCA(m1_rateM))
## $Subj
## Importance of components:
##                          [,1]   [,2]    [,3] [,4]
## Standard deviation     3.5108 2.5318 0.55486    0
## Proportion of Variance 0.6472 0.3366 0.01617    0
## Cumulative Proportion  0.6472 0.9838 1.00000    1
## 
## $item
## Importance of components:
##                          [,1]
## Standard deviation     0.1407
## Proportion of Variance 1.0000
## Cumulative Proportion  1.0000
# 1 of 4 components for Subj has zero variance | intercept only for item
### ANOVA model comparison - initial model - reduced model #############################
anova(m0_rateM, m1_rateM)
# not significant, AIC smaller for m1_rateM - preferred
### SECOND #############################
# try exclusion of randEff for item altogether
summary(m2_rateM <- glmer(rate ~ 1+
                            grpd.ungr+lhd.rhd+rd.rpt+
                            grpd.ungr_lhd.rhd+grpd.ungr_rd.rpt+lhd.rhd_rd.rpt+
                            grpd.ungr_lhd.rhd_rd.rpt+
                            (1+grpd.ungr+rd.rpt+grpd.ungr_rd.rpt|Subj),
                          family="binomial",
                          control=glmerControl(optCtrl=list(maxfun=1e5),
                                               optimizer="bobyqa"),
                          pDat_mm))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: rate ~ 1 + grpd.ungr + lhd.rhd + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_rd.rpt + lhd.rhd_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     (1 + grpd.ungr + rd.rpt + grpd.ungr_rd.rpt | Subj)
##    Data: pDat_mm
## Control: glmerControl(optCtrl = list(maxfun = 1e+05), optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1168.4   1261.3   -566.2   1132.4     1271 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7248 -0.1564  0.2837  0.4990  3.1065 
## 
## Random effects:
##  Groups Name             Variance Std.Dev. Corr          
##  Subj   (Intercept)      2.523    1.588                  
##         grpd.ungr        8.062    2.839    0.95          
##         rd.rpt           4.045    2.011    0.34 0.38     
##         grpd.ungr_rd.rpt 4.295    2.073    0.12 0.24 0.95
## Number of obs: 1289, groups:  Subj, 32
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.5686     0.3588   4.371 1.24e-05 ***
## grpd.ungr                  0.1852     0.6553   0.283 0.777463    
## lhd.rhd                    0.7943     0.6614   1.201 0.229802    
## rd.rpt                     2.3150     0.5632   4.110 3.95e-05 ***
## grpd.ungr_lhd.rhd          1.0815     1.1858   0.912 0.361757    
## grpd.ungr_rd.rpt           3.4936     0.9221   3.789 0.000151 ***
## lhd.rhd_rd.rpt             1.4963     0.9872   1.516 0.129608    
## grpd.ungr_lhd.rhd_rd.rpt   2.2828     1.4906   1.531 0.125661    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grpd.n lhd.rh rd.rpt grpd.ngr_l. grpd.ngr_r. lhd._.
## grpd.ungr   0.908                                                     
## lhd.rhd     0.261  0.242                                              
## rd.rpt      0.433  0.461  0.148                                       
## grpd.ngr_l. 0.245  0.271  0.901  0.170                                
## grpd.ngr_r. 0.363  0.427  0.183  0.832  0.202                         
## lhd.rhd_rd. 0.153  0.166  0.291  0.374  0.330       0.375             
## grpd.ng_._. 0.195  0.207  0.177  0.407  0.236       0.514       0.796 
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
summary(rePCA(m2_rateM))
## $Subj
## Importance of components:
##                          [,1]   [,2]    [,3] [,4]
## Standard deviation     3.4996 2.5243 0.55304    0
## Proportion of Variance 0.6471 0.3367 0.01616    0
## Cumulative Proportion  0.6471 0.9838 1.00000    1
# 1 of 4 components for Subj has zero variance
### ANOVA model comparison - m2_rateM - m1_rateM #############################
anova(m1_rateM, m2_rateM)
# not significant, AIC smaller for m2_rateM - preferred
### THIRD #############################
# try exclusion of "rd.rpt" for Subj
summary(m3_rateM <- glmer(rate ~ 1+
                            grpd.ungr+lhd.rhd+rd.rpt+
                            grpd.ungr_lhd.rhd+grpd.ungr_rd.rpt+lhd.rhd_rd.rpt+
                            grpd.ungr_lhd.rhd_rd.rpt+
                            (1+grpd.ungr+grpd.ungr_rd.rpt|Subj),
                          family="binomial",
                          control=glmerControl(optCtrl=list(maxfun=1e5),
                                               optimizer="bobyqa"),
                          pDat_mm))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: rate ~ 1 + grpd.ungr + lhd.rhd + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_rd.rpt + lhd.rhd_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     (1 + grpd.ungr + grpd.ungr_rd.rpt | Subj)
##    Data: pDat_mm
## Control: glmerControl(optCtrl = list(maxfun = 1e+05), optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1204.2   1276.5   -588.1   1176.2     1275 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.9531 -0.2444  0.3198  0.5289  2.9811 
## 
## Random effects:
##  Groups Name             Variance Std.Dev. Corr       
##  Subj   (Intercept)      1.962    1.401               
##         grpd.ungr        5.799    2.408     0.91      
##         grpd.ungr_rd.rpt 3.459    1.860    -0.53 -0.22
## Number of obs: 1289, groups:  Subj, 32
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.3083     0.2767   4.728 2.27e-06 ***
## grpd.ungr                 -0.3525     0.4953  -0.712 0.476644    
## lhd.rhd                    0.4416     0.5474   0.807 0.419790    
## rd.rpt                     1.8220     0.2307   7.899 2.82e-15 ***
## grpd.ungr_lhd.rhd          0.5074     0.9728   0.522 0.601954    
## grpd.ungr_rd.rpt           2.2157     0.5788   3.828 0.000129 ***
## lhd.rhd_rd.rpt             1.3443     0.4454   3.018 0.002543 ** 
## grpd.ungr_lhd.rhd_rd.rpt   1.9777     1.1119   1.779 0.075286 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grpd.n lhd.rh rd.rpt grpd.ngr_l. grpd.ngr_r. lhd._.
## grpd.ungr    0.822                                                    
## lhd.rhd      0.167  0.148                                             
## rd.rpt       0.069  0.025  0.019                                      
## grpd.ngr_l.  0.145  0.193  0.835  0.026                               
## grpd.ngr_r. -0.239 -0.033  0.007  0.370  0.022                        
## lhd.rhd_rd.  0.045  0.046  0.028  0.542  0.012       0.336            
## grpd.ng_._.  0.004  0.022 -0.279  0.317 -0.097       0.394       0.399
summary(rePCA(m3_rateM))
## $Subj
## Importance of components:
##                          [,1]   [,2]    [,3]
## Standard deviation     2.8351 1.7611 0.28357
## Proportion of Variance 0.7164 0.2764 0.00717
## Cumulative Proportion  0.7164 0.9928 1.00000
# all components seem okay
### ANOVA model comparison - m3_rateM - m2_rateM #############################
anova(m2_rateM, m3_rateM)
# significant, AIC smaller for m2_rateM - preferred
### FOURTH #############################
# back to m2_rateM
# try exclusion of "grpd.ungr_rd.rpt" for Subj
summary(m4_rateM <- glmer(rate ~ 1+
                            grpd.ungr+lhd.rhd+rd.rpt+
                            grpd.ungr_lhd.rhd+grpd.ungr_rd.rpt+lhd.rhd_rd.rpt+
                            grpd.ungr_lhd.rhd_rd.rpt+
                            (1+grpd.ungr+rd.rpt|Subj),
                          family="binomial",
                          control=glmerControl(optCtrl=list(maxfun=1e5),
                                               optimizer="bobyqa"),
                          pDat_mm))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: rate ~ 1 + grpd.ungr + lhd.rhd + rd.rpt + grpd.ungr_lhd.rhd +  
##     grpd.ungr_rd.rpt + lhd.rhd_rd.rpt + grpd.ungr_lhd.rhd_rd.rpt +  
##     (1 + grpd.ungr + rd.rpt | Subj)
##    Data: pDat_mm
## Control: glmerControl(optCtrl = list(maxfun = 1e+05), optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1175.8   1248.1   -573.9   1147.8     1275 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9809 -0.1912  0.2966  0.4746  3.0101 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr     
##  Subj   (Intercept) 2.098    1.449             
##         grpd.ungr   6.509    2.551    0.93     
##         rd.rpt      2.417    1.555    0.27 0.34
## Number of obs: 1289, groups:  Subj, 32
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.4875     0.2950   5.042 4.60e-07 ***
## grpd.ungr                 -0.1468     0.5171  -0.284    0.777    
## lhd.rhd                    0.6953     0.5797   1.199    0.230    
## rd.rpt                     1.9285     0.3983   4.842 1.29e-06 ***
## grpd.ungr_lhd.rhd          0.8638     1.0209   0.846    0.397    
## grpd.ungr_rd.rpt           2.9187     0.5246   5.563 2.65e-08 ***
## lhd.rhd_rd.rpt             1.3290     0.7716   1.722    0.085 .  
## grpd.ungr_lhd.rhd_rd.rpt   2.0324     0.9946   2.043    0.041 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grpd.n lhd.rh rd.rpt grpd.ngr_l. grpd.ngr_r. lhd._.
## grpd.ungr   0.835                                                     
## lhd.rhd     0.188  0.146                                              
## rd.rpt      0.247  0.291  0.051                                       
## grpd.ngr_l. 0.147  0.181  0.842  0.083                                
## grpd.ngr_r. 0.138  0.108  0.078  0.382  0.071                         
## lhd.rhd_rd. 0.065  0.088  0.194  0.350  0.274       0.306             
## grpd.ng_._. 0.107  0.101  0.092  0.326  0.074       0.600       0.374
summary(rePCA(m4_rateM))
## $Subj
## Importance of components:
##                          [,1]   [,2]   [,3]
## Standard deviation     2.9565 1.4417 0.4529
## Proportion of Variance 0.7929 0.1885 0.0186
## Cumulative Proportion  0.7929 0.9814 1.0000
# all components seem okay
### ANOVA model comparison - m4_rateM - m2_rateM #############################
anova(m4_rateM, m2_rateM)
# significant, AIC smaller for m2_rateM - preferred
# no better model fit achievable
# use m2_rateM
# FINAL model ---------------------------------
## model with factors and * for fixed effects =================================
# use REML=TRUE
summary(rateM <- glmer(rate ~ condition*group*exp+
                         (1+grpd.ungr+rd.rpt+grpd.ungr_rd.rpt|Subj),
                       family="binomial",
                       control=glmerControl(optCtrl=list(maxfun=1e5),
                                            optimizer="bobyqa"),
                       pDat_mm))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## rate ~ condition * group * exp + (1 + grpd.ungr + rd.rpt + grpd.ungr_rd.rpt |  
##     Subj)
##    Data: pDat_mm
## Control: glmerControl(optCtrl = list(maxfun = 1e+05), optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1168.4   1261.3   -566.2   1132.4     1271 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7248 -0.1564  0.2837  0.4990  3.1065 
## 
## Random effects:
##  Groups Name             Variance Std.Dev. Corr          
##  Subj   (Intercept)      2.523    1.588                  
##         grpd.ungr        8.062    2.839    0.95          
##         rd.rpt           4.045    2.011    0.34 0.38     
##         grpd.ungr_rd.rpt 4.295    2.073    0.12 0.24 0.95
## Number of obs: 1289, groups:  Subj, 32
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              1.5686     0.3588   4.371 1.24e-05 ***
## condition1               0.1852     0.6553   0.283 0.777463    
## group1                   0.7943     0.6614   1.201 0.229802    
## exp1                     2.3150     0.5632   4.110 3.95e-05 ***
## condition1:group1        1.0815     1.1858   0.912 0.361757    
## condition1:exp1          3.4936     0.9221   3.789 0.000151 ***
## group1:exp1              1.4963     0.9872   1.516 0.129608    
## condition1:group1:exp1   2.2828     1.4906   1.531 0.125661    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtn1 group1 exp1  cndtn1:g1 cndtn1:x1 grp1:1
## condition1  0.908                                                
## group1      0.261  0.242                                         
## exp1        0.433  0.461  0.148                                  
## cndtn1:grp1 0.245  0.271  0.901  0.170                           
## condtn1:xp1 0.363  0.427  0.183  0.832 0.202                     
## group1:exp1 0.153  0.166  0.291  0.374 0.330     0.375           
## cndtn1:g1:1 0.195  0.207  0.177  0.407 0.236     0.514     0.796 
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# isSingular
save(rateM, file = "rateM.Rdata")
# *Nice table output* ---------------------------------
# sjPlot::tab_model((rateM), show.se = TRUE, show.stat = TRUE, collapse.ci=TRUE, transform = NULL, file="rateM.doc")
# clean workspace ---------------------------------
# ls()
rm(m0_rateM, m1_rateM, m2_rateM, m3_rateM, m4_rateM, zcp_rateM)

9 Production accuracy model

9.1 Visualizations and Post-Hoc-Tests

# interaction visualization ---------------------------------
load("rateM.Rdata")
# condition:exp =================================
### condition ~ exp #############################
emmeans::emmip(rateM, condition ~ exp, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - production accuracy",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

### exp ~ condition #############################
emmeans::emmip(rateM, exp ~ condition, CIs = TRUE, 
               CIarg = list(lwd = 1.1, alpha=0.7), 
               linearg = list(linetype = "dashed", alpha=0)) +
  labs(
    y = "Model predictions - production accuracy",
  ) +
  scale_color_brewer(palette = "Set1") +
  theme(axis.text=element_text(size=11),
        axis.title.x=element_text(size=11),
        axis.title.y=element_text(size=11),
        legend.justification = c("right", "top"),
        legend.title = element_text(size=11, face="bold"), 
        legend.text=element_text(size=11, face="bold")
  )

# POST-HOC testing ---------------------------------
# what are the significant differences
emmeans::emmeans(rateM, pairwise ~ condition*exp) %>% 
    summary(by = NULL, adjust = "bonferroni")
## $emmeans
##  condition exp           emmean    SE  df asymp.LCL asymp.UCL
##  grouped   reading_aloud   3.69 0.994 Inf     1.210      6.17
##  ungrouped reading_aloud   1.76 0.232 Inf     1.181      2.34
##  grouped   repetition     -0.37 0.627 Inf    -1.936      1.20
##  ungrouped repetition      1.19 0.200 Inf     0.693      1.69
## 
## Results are averaged over the levels of: group 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## 
## $contrasts
##  contrast                                        estimate    SE  df z.ratio
##  grouped reading_aloud - ungrouped reading_aloud    1.932 0.949 Inf   2.037
##  grouped reading_aloud - grouped repetition         4.062 0.981 Inf   4.141
##  grouped reading_aloud - ungrouped repetition       2.500 1.043 Inf   2.398
##  ungrouped reading_aloud - grouped repetition       2.130 0.637 Inf   3.342
##  ungrouped reading_aloud - ungrouped repetition     0.568 0.312 Inf   1.820
##  grouped repetition - ungrouped repetition         -1.562 0.620 Inf  -2.520
##  p.value
##   0.2501
##   0.0002
##   0.0989
##   0.0050
##   0.4124
##   0.0705
## 
## Results are averaged over the levels of: group 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: bonferroni method for 6 tests

10 Raw data plots

10.1 Production accuracy

10.2 Plotting means and error bars

Functions for mean + CI plots
+ start with this script chunk “SEsummary” + run it and then go on to the next chunk to plot the data on all cues “cues separately”

# summarySE between ---------------------------------
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%)
## data:          a data frame
## measurevar:    the name of a column that contains the variable to be summarized
## groupvars:     a vector containing names of columns that contain grouping variables
## na.rm:         a boolean that indicates whether to ignore NA's
## conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=TRUE,
                      conf.interval=.95, .drop=TRUE) {
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=TRUE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  # Rename the "mean" column    
  datac <- plyr::rename(datac, c("mean" = measurevar))
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  return(datac)
}
# normDataWithin ---------------------------------
## Norms the data within specified groups in a data frame; it normalizes each
## subject (identified by idvar) so that they have the same mean, within each group
## specified by betweenvars
## data:        a data frame.
## idvar:       the name of a column that identifies each subject (or matched subjects)
## measurevar:  the name of a column that contains the variable to be summariezed
## betweenvars: a vector containing names of columns that are between-subjects variables
## na.rm:       a boolean that indicates whether to ignore NA's
normDataWithin <- function(data=NULL, idvar, measurevar, betweenvars=NULL,
                           na.rm=TRUE, .drop=TRUE) {
  
  # Measure var on left, idvar + between vars on right of formula
  data.subjMean <- ddply(data, c(idvar, betweenvars), .drop=.drop,
                         .fun = function(xx, col, na.rm) {
                           c(subjMean = mean(xx[,col], na.rm=na.rm))
                         },
                         measurevar,
                         na.rm
  )
  
  # Put the subject means with original data
  data <- merge(data, data.subjMean)
  
  # Get the normalized data in a new column
  measureNormedVar <- paste(measurevar, "_norm", sep="")
  data[,measureNormedVar] <- data[,measurevar] - data[,"subjMean"] +
    mean(data[,measurevar], na.rm=na.rm)
  
  # Remove this subject mean column
  data$subjMean <- NULL
  
  return(data)
}
# SEsummarywithin ---------------------------------
## Summarizes data, handling within-subjects variables by removing inter-subject variability.
## It will still work if there are no within-S variables.
## Gives count, un-normed mean, normed mean (with same between-group mean),
## standard deviation, standard error of the mean, and confidence interval.
## If there are within-subject variables, calculate adjusted values using method from Morey (2008)
## data:          a data frame.
## measurevar:    the name of a column that contains the variable to be summariezed
## betweenvars:   a vector containing names of columns that are between-subjects variables
## withinvars:    a vector containing names of columns that are within-subjects variables
## idvar:         the name of a column that identifies each subject (or matched subjects)
## na.rm:         a boolean that indicates whether to ignore NA's
## conf.interval: the percent range of the confidence interval (default is 95%)
summarySEwithin <- function(data=NULL, measurevar, betweenvars=NULL, withinvars=NULL,
                            idvar=NULL, na.rm=TRUE, conf.interval=.95, .drop=TRUE) {
  
  # Ensure that the betweenvars and withinvars are factors
  factorvars <- vapply(data[, c(betweenvars, withinvars), drop=FALSE],
                       FUN=is.factor, FUN.VALUE=logical(1))
  
  if (!all(factorvars)) {
    nonfactorvars <- names(factorvars)[!factorvars]
    message("Automatically converting the following non-factors to factors: ",
            paste(nonfactorvars, collapse = ", "))
    data[nonfactorvars] <- lapply(data[nonfactorvars], factor)
  }
  
  # Get the means from the un-normed data
  datac <- summarySE(data, measurevar, groupvars=c(betweenvars, withinvars),
                     na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)
  
  # Drop all the unused columns (these will be calculated with normed data)
  datac$sd <- NULL
  datac$se <- NULL
  datac$ci <- NULL
  
  # Norm each subject's data
  ndata <- normDataWithin(data, idvar, measurevar, betweenvars, na.rm, .drop=.drop)
  
  # This is the name of the new column
  measurevar_n <- paste(measurevar, "_norm", sep="")
  
  # Collapse the normed data - now we can treat between and within vars the same
  ndatac <- summarySE(ndata, measurevar_n, groupvars=c(betweenvars, withinvars),
                      na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)
  
  # Apply correction from Morey (2008) to the standard error and confidence interval
  #  Get the product of the number of conditions of within-S variables
  nWithinGroups    <- prod(vapply(ndatac[,withinvars, drop=FALSE], FUN=nlevels,
                                  FUN.VALUE=numeric(1)))
  correctionFactor <- sqrt( nWithinGroups / (nWithinGroups-1) )
  
  # Apply the correction factor
  ndatac$sd <- ndatac$sd * correctionFactor
  ndatac$se <- ndatac$se * correctionFactor
  ndatac$ci <- ndatac$ci * correctionFactor
  
  # Combine the un-normed means with the normed results
  merge(datac, ndatac)
}

Mean an CI plot for cues on name2

# Cues ---------------------------------
Cues <- pDat %>% 
  dplyr::select(Subj, item, group, exp, condition, accuracy, f0_range2, s8reln2, pause3) %>% 
  tidyr::gather("variable", "measurement", -Subj, -item, -group, -exp, -condition, -accuracy)
Cues$group <- as.factor(as.character(Cues$group))
Cues$exp <- as.factor(as.character(Cues$exp))
Cues$condition <- as.factor(as.character(Cues$condition))
Cues$accuracy <- as.factor(as.character(Cues$accuracy))
Cues$variable <- as.factor(as.character(Cues$variable))
Cues$measurement <- as.numeric(as.character(Cues$measurement))
## f0_range =================================
f0.range <- Cues %>% 
  filter(variable=="f0_range2") %>% 
  group_by(Subj, group, exp, condition, accuracy, variable) %>% 
  summarySEwithin(measurevar = "measurement",
                  betweenvars = "group",
                  withinvars = c("variable", "condition", "exp", "accuracy"),
                  idvar = "Subj") %>% 
  ggplot(aes(x=accuracy,y=measurement, shape = condition, color=condition)) +
  geom_point(position = position_dodge2(width = .5)) + 
  facet_grid(exp~group)+
  geom_errorbar(aes(ymax = measurement + ci, ymin = measurement - ci), position = position_dodge2(width = .2), width=.5) +
  theme_bw() +
  theme(axis.text=element_text(size=14), 
        axis.title=element_text(size=14),
        strip.text.x = element_text(size=14),
        strip.text.y = element_blank(),
        legend.title=element_text(size=14), 
        legend.text=element_text(size=14)) +
  labs(
    x = "",
    y = "f0-range on name2 (st)"
  )
f0.range

# ggsave("f0_range2.jpg", width = 18, height = 14, dpi = 300, units = "cm")
## final lengthening =================================
fin.length <- Cues %>% 
  filter(variable=="s8reln2") %>% 
  group_by(Subj, group, exp, condition, accuracy, variable) %>% 
  summarySEwithin(measurevar = "measurement",
                  betweenvars = "group",
                  withinvars = c("variable", "condition", "exp", "accuracy"),
                  idvar = "Subj") %>% 
  ggplot(aes(x=accuracy,y=measurement, shape = condition, color=condition)) +
  geom_point(position = position_dodge2(width = .5)) + 
  facet_grid(exp~group)+
  geom_errorbar(aes(ymax = measurement + ci, ymin = measurement - ci), position = position_dodge2(width = .2), width=.5) +
  theme_bw() +
  theme(axis.text=element_text(size=14), 
        axis.title=element_text(size=14),
        strip.text.x = element_text(size=14),
        strip.text.y = element_blank(),
        legend.title=element_text(size=14), 
        legend.text=element_text(size=14)) +
  labs(
    x = "",
    y = "final lengthening on name2 (%)"
  )
fin.length

# ggsave("fin-length.jpg", width = 18, height = 14, dpi = 300, units = "cm")
## pause 3 =================================
pause3 <- Cues %>% 
  filter(variable=="pause3") %>% 
  group_by(Subj, group, exp, condition, accuracy, variable) %>% 
  summarySEwithin(measurevar = "measurement",
                  betweenvars = "group",
                  withinvars = c("variable", "condition", "exp", "accuracy"),
                  idvar = "Subj") %>% 
  ggplot(aes(x=accuracy,y=measurement, shape = condition, color=condition)) +
  geom_point(position = position_dodge2(width = .5)) + 
  facet_grid(exp~group)+
  geom_errorbar(aes(ymax = measurement + ci, ymin = measurement - ci), position = position_dodge2(width = .2), width=.5) +
  theme_bw() +
  theme(axis.text=element_text(size=14), 
        axis.title=element_text(size=14),
        strip.text.x = element_text(size=14),
        strip.text.y = element_text(size=14),
        legend.title=element_text(size=14), 
        legend.text=element_text(size=14)) +
  labs(
    x = "",
    y = "pause after name2 (%)"
  )
pause3

# ggsave("pause3.jpg", width = 18, height = 14, dpi = 300, units = "cm")
## plot grids =================================
library(ggpubr)
pDat_cues <- ggarrange(f0.range,fin.length,pause3, ncol=3, nrow=1, common.legend = TRUE, legend="bottom")
pDat_cues

# ggsave("cue_grid.jpg", plot = pDat_cues, width = 18, height = 10, dpi = 300)
## clean workspace =================================
rm(Cues, f0.range, fin.length, pause3, pDat_cues, normDataWithin, summarySE, summarySEwithin)

11 Correlations - working memory measures / accuracy

# import working memory data ---------------------------------
workMem <- read.csv2("data/workMem.csv", header = TRUE)
# enumerate accuracy variable ---------------------------------
pDat$rating_match <- as.numeric(as.factor(pDat$rating_match)) -1
## participant check =================================
## list first entry for exp for each participant
pDat %>%
  group_by(Subj, group, exp) %>%
  dplyr::summarise(mean = round(mean(rating_match)/11*100, digits=2)) %>% 
  group_by(Subj) %>%
  summarise(exp = first(exp)) %>%
  as.data.frame()
### RHDP11 only took part in repetition exp
## list how many entries for exp for each participant exist
pDat %>%
  group_by(Subj, group, exp) %>%
  dplyr::summarise(mean = round(mean(rating_match)/11*100, digits=2)) %>% 
  group_by(Subj) %>%
  summarise(exp = n_distinct(exp)) %>%
  as.data.frame()
# LHDP01, LHDP02, LHDP08, LHDP12, RHDP19, RHDP20 only took part in reading_aloud
## working memory subset =================================
## only individuals that participated in both exp
workMemBoth <- subset(workMem, 
                      workMem$Subj=="LHDP03" | workMem$Subj=="LHDP04" | workMem$Subj=="LHDP07" | workMem$Subj=="LHDP09" |
                        workMem$Subj=="LHDP10" | workMem$Subj=="LHDP11" | workMem$Subj=="LHDP13" | workMem$Subj=="LHDP14" |
                        workMem$Subj=="LHDP15" | workMem$Subj=="LHDP16" | workMem$Subj=="RHDP01" | workMem$Subj=="RHDP02" |
                        workMem$Subj=="RHDP03" | workMem$Subj=="RHDP04" | workMem$Subj=="RHDP05" | workMem$Subj=="RHDP08" |
                        workMem$Subj=="RHDP09" | workMem$Subj=="RHDP12" | workMem$Subj=="RHDP13" | workMem$Subj=="RHDP14" |
                        workMem$Subj=="RHDP15" | workMem$Subj=="RHDP16" | workMem$Subj=="RHDP17" | workMem$Subj=="RHDP18" |
                        workMem$Subj=="RHDP22")
## data subset =================================
## only individuals that participated in both exp
meanBoth <-  pDat %>% 
  group_by(Subj) %>% 
  filter(n_distinct(exp)==2)
# data aggregation  ---------------------------------
### means per Subj, group, exp, condition merged with working memory measures  #############################
mean_rating_Subj <- meanBoth %>%
  group_by(Subj, group, exp, condition) %>%
  dplyr::summarise(mean = round(mean(rating_match)/11*100, digits=2))
allMeans <- merge(mean_rating_Subj, workMemBoth[, c("Subj", "compDigit")], by=c("Subj"), all.x=TRUE)
### means per Subj, group, exp - overall conditions merged with working memory measures  #############################
mean_rating_Subj1 <- meanBoth %>%
  group_by(Subj, group, exp) %>%
  dplyr::summarise(mean = round(mean(rating_match)/11*100, digits=2))
allMeans1 <- merge(mean_rating_Subj1, workMemBoth[, c("Subj", "compDigit")], by=c("Subj"), all.x=TRUE)
# plots ---------------------------------
### compDigit - COMPOSITE SCORE DIGIT SPAN forward & backward #############################
## REPETITION =================================
# mean accuracy assessment over group and OVERALL conditions
### these correlational measures are REPORTED IN PAPER ("3.2.3 Effects of Working Memory") ###
corrRep <- facet(ggscatter(subset(allMeans1, allMeans1$exp=="repetition"), x="mean", y="compDigit", 
                           add = "reg.line", conf.int = TRUE, 
                           cor.coef = TRUE, cor.coeff.args = list(size=6, label.x = 50, label.y = 80, label.sep = "\n"), cor.method = "spearman",
                           xlab = "mean accuracy assessment repetition overall", ylab = "composite score digit span")+
                   theme_bw(),
                 facet.by = c("group"),
                 short.panel.labs = TRUE, # Allow long labels in panels
                 panel.labs.background = list(fill = "lightgrey", color = "lightgrey"),
                 panel.labs.font = list(face = NULL, color = "black", size = 20, angle = NULL))+
  font("xlab", size = 20, color = "black")+
  font("ylab", size = 20, color = "black")+
  font("xy.text", size = 19, color = "black")
corrRep

# ggsave("correlations_repetition_overall.jpg", plot = corrRep1, width = 18, height = 16, dpi = 300)
# mean accuracy assessment over condition AND group
corrRep1 <- facet(ggscatter(subset(allMeans, allMeans$exp=="repetition"), x="mean", y="compDigit", 
                            add = "reg.line", conf.int = TRUE, 
                            cor.coef = TRUE, cor.coeff.args = list(size=6, label.x = 25, label.y = 75, label.sep = "\n"), cor.method = "spearman",
                            xlab = "mean accuracy assessment repetition", ylab = "composite score digit span")+
                    theme_bw(),
                  facet.by = c("group","condition"),
                  short.panel.labs = TRUE, # Allow long labels in panels
                  panel.labs.background = list(fill = "lightgrey", color = "lightgrey"),
                  panel.labs.font = list(face = NULL, color = "black", size = 20, angle = NULL))+
  font("xlab", size = 20, color = "black")+
  font("ylab", size = 20, color = "black")+
  font("xy.text", size = 19, color = "black")
corrRep1

# ggsave("correlations_repetition.jpg", plot = corrRep, width = 18, height = 16, dpi = 300)
# READING ALOUD  =================================
# mean accuracy assessment over group and OVERALL conditions
corrRead <- facet(ggscatter(subset(allMeans1, allMeans1$exp=="reading_aloud"), x="mean", y="compDigit", 
                            add = "reg.line", conf.int = TRUE, 
                            cor.coef = TRUE, cor.coeff.args = list(size=6, label.x = 50, label.y = 80, label.sep = "\n"), cor.method = "spearman",
                            xlab = "mean accuracy assessment reading_aloud overall", ylab = "composite score digit span")+
                    theme_bw(),
                  facet.by = c("group"),
                  short.panel.labs = TRUE, # Allow long labels in panels
                  panel.labs.background = list(fill = "lightgrey", color = "lightgrey"),
                  panel.labs.font = list(face = NULL, color = "black", size = 20, angle = NULL))+
  font("xlab", size = 20, color = "black")+
  font("ylab", size = 20, color = "black")+
  font("xy.text", size = 19, color = "black")
corrRead

# ggsave("correlations_reading_aloud_overall.jpg", plot = corrRead1, width = 18, height = 16, dpi = 300)
# mean accuracy assessment over condition AND group
corrRead1 <- facet(ggscatter(subset(allMeans, allMeans$exp=="reading_aloud"), x="mean", y="compDigit", 
                             add = "reg.line", conf.int = TRUE, 
                             cor.coef = TRUE, cor.coeff.args = list(size=6, label.x = 20, label.y = 75, label.sep = "\n"), cor.method = "spearman",
                             xlab = "mean accuracy assessment reading_aloud", ylab = "composite score digit span")+
                     theme_bw(),
                   facet.by = c("group","condition"),
                   short.panel.labs = TRUE, # Allow long labels in panels
                   panel.labs.background = list(fill = "lightgrey", color = "lightgrey"),
                   panel.labs.font = list(face = NULL, color = "black", size = 20, angle = NULL))+
  font("xlab", size = 20, color = "black")+
  font("ylab", size = 20, color = "black")+
  font("xy.text", size = 19, color = "black")
corrRead1

# ggsave("correlations_reading_aloud.jpg", plot = corrRead, width = 18, height = 16, dpi = 300)

12 Summary Tables of Cue Use for Patient Groups and Model-Stimuli

require(flextable)
data_corr <- subset(data, data$accuracy=="correct")
# generate descriptive summary for pause 3 per group (pause3 duration relative to utterance duration in %) ---------------------------------
pause3group <- data_corr %>%
  dplyr::group_by(exp, group, condition) %>%
  dplyr::summarise(mean = round(mean(pause3),3), n = n(), min = round(min(pause3),3), 
                   max = round(max(pause3),3), sd = round(sd(pause3),3))
## generate and safe respective table for pause 3 per group =================================
pause3groupFT <- flextable(pause3group) %>%
  merge_v(j = c("exp", "group")) %>% 
  valign(valign = "top") %>% 
  theme_box() %>% set_table_properties(layout = "autofit")
pause3groupFT <- add_header_row(
  x = pause3groupFT, values = c("Descriptive summary for pause3 by patient group and model stimuli
                                   Only correct production accuracy"),
  colwidths = c(8))
pause3groupFT <- align(pause3groupFT, i = 1, part = "header", align = "center")
pause3groupFT <- bold(pause3groupFT, j = 1, i = ~ !is.na(exp))
pause3groupFT <- fontsize(pause3groupFT, part = "all", size = 12)
pause3groupFT <- theme_vanilla(pause3groupFT)
pause3groupFT
save_as_docx(
  "Table Summary Pause Cue" = pause3groupFT, 
  path = "group_sum_pause.docx")
rm(pause3group)
rm(pause3groupFT)
# F0-range ---------------------------------
# generate descriptive summary for f0_range2 per group (difference in semitones btw L and H annotated on name2) ---------------------------------
f0_range2group <- data_corr %>%
  dplyr::group_by(exp, group, condition) %>%
  dplyr::summarise(mean = round(mean(f0_range2, na.rm=T),3), n = n(), miss = sum(is.na(f0_range2)), min = round(min(f0_range2, na.rm=T),3), 
                   max = round(max(f0_range2, na.rm=T),3), sd = round(sd(f0_range2, na.rm=T),3))
## generate and safe respective table for f0_range2 per group =================================
f0_range2groupFT <- flextable(f0_range2group) %>%
  merge_v(j = c("exp", "group")) %>% 
  valign(valign = "top") %>% 
  theme_box() %>% set_table_properties(layout = "autofit")
f0_range2groupFT <- add_header_row(
  x = f0_range2groupFT, values = c("Descriptive summary for f0-range on name2 by patient group and model stimuli
                                   Only correct production accuracy"),
  colwidths = c(9))
f0_range2groupFT <- align(f0_range2groupFT, i = 1, part = "header", align = "center")
f0_range2groupFT <- bold(f0_range2groupFT, j = 1, i = ~ !is.na(exp))
f0_range2groupFT <- fontsize(f0_range2groupFT, part = "all", size = 12)
f0_range2groupFT <- theme_vanilla(f0_range2groupFT)
f0_range2groupFT
save_as_docx(
  "Table Summary f0-Range Cue" = f0_range2groupFT, 
  path = "group_sum_f0.docx")
rm(f0_range2group)
rm(f0_range2groupFT)
# Final lengthening ---------------------------------
# generate descriptive summary for s8reln2 (length of final vowel, s8, in % relative to duration of name2) ---------------------------------
s8reln2group <- data_corr %>%
  dplyr::group_by(exp, group, condition) %>%
  dplyr::summarise(mean = round(mean(s8reln2),3), n = n(), min = round(min(s8reln2),3), 
                   max = round(max(s8reln2),3), sd = round(sd(s8reln2),3))
## 8B generate and safe respective table for s8reln2  =================================
s8reln2groupFT <- flextable(s8reln2group) %>%
  merge_v(j = c("exp", "group")) %>% 
  valign(valign = "top") %>% 
  theme_box() %>% set_table_properties(layout = "autofit")
s8reln2groupFT <- add_header_row(
  x = s8reln2groupFT, values = c("Descriptive summary for final lengthening on name2 by patient group and model stimuli
                                   Only correct production accuracy"),
  colwidths = c(8))
s8reln2groupFT <- align(s8reln2groupFT, i = 1, part = "header", align = "center")
s8reln2groupFT <- bold(s8reln2groupFT, j = 1, i = ~ !is.na(exp))
s8reln2groupFT <- fontsize(s8reln2groupFT, part = "all", size = 12)
s8reln2groupFT <- theme_vanilla(s8reln2groupFT)
s8reln2groupFT
save_as_docx(
  "Table Summary Final Lengthening Cue" = s8reln2groupFT, 
  path = "group_sum_FL.docx")
rm(s8reln2group)
rm(s8reln2groupFT)

13 Boxplot model-stimuli

mod_stim <- subset(data_corr, data_corr$Subj=="model_stimulus")
# cue value distribution by condition significantly different from each other?
mod_grou <- subset(mod_stim, mod_stim$condition=="grouped")
mod_ungrou <- subset(mod_stim, mod_stim$condition=="ungrouped")
# vectors / groups for cues
x1 <- mod_grou$f0_range2
x2 <- mod_ungrou$f0_range2

y1 <- mod_grou$s8reln2
y2 <- mod_ungrou$s8reln2

z1 <- mod_grou$pause3
z2 <- mod_ungrou$pause3
# dependent 2-group Wilcoxon Signed Rank Test
# where both groups are numeric
wilcox.test(x1,x2,paired=TRUE) 
## 
##  Wilcoxon signed rank exact test
## 
## data:  x1 and x2
## V = 21, p-value = 0.03125
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(y1,y2,paired=TRUE)
## 
##  Wilcoxon signed rank exact test
## 
## data:  y1 and y2
## V = 21, p-value = 0.03125
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(z1,z2,paired=TRUE)
## 
##  Wilcoxon signed rank exact test
## 
## data:  z1 and z2
## V = 21, p-value = 0.03125
## alternative hypothesis: true location shift is not equal to 0
# boxplot f0-range2
f0 <- ggplot(mod_stim, aes(condition, f0_range2))
f0_p <- f0 + geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.1) +
  labs(
    title = "F0 range on name2 by condition for model stimuli",
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14)) +
  annotate("segment", x = 1, xend = 2, y = 11.5, yend = 11.5,
  colour = "black") +
  annotate("text", x = 1.5, y = 11, label = "Wilcoxon test
           p-value = .031 (V=21)")
# boxplot s8reln2
fl <- ggplot(mod_stim, aes(condition, s8reln2))
fl_p <- fl + geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.1) +
  labs(
    title = "Final Lengthening on name2 by condition for model stimuli",
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14)) +
  annotate("segment", x = 1, xend = 2, y = 49, yend = 49,
  colour = "black") +
  annotate("text", x = 1.5, y = 48, label = "Wilcoxon test
           p-value = .031 (V=21)")
# boxplot pause3
p3 <- ggplot(mod_stim, aes(condition, pause3))
p3_p <- p3 + geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.1)+
  labs(
    title = "Pause after name2 by condition for model stimuli",
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14)) +
  annotate("segment", x = 1, xend = 2, y = 24, yend = 24,
  colour = "black") +
  annotate("text", x = 1.5, y = 23, label = "Wilcoxon test
           p-value = .031 (V=21)")
## plot grids =================================
require(ggpubr)
mod_stim_cues <- ggarrange(f0_p,fl_p,p3_p, ncol=3, nrow=1)
mod_stim_cues

# ggsave("mod_stim_cue.jpg", plot = mod_stim_cues, width = 18, height = 10, dpi = 300)

14 Session information

sessionInfo(package = NULL)
## R version 4.2.1 (2022-06-23)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.5
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] flextable_0.7.2 ggpubr_0.4.0    lmerTest_3.1-3  lme4_1.1-30    
##  [5] Matrix_1.4-1    forcats_0.5.1   stringr_1.4.0   dplyr_1.0.9    
##  [9] purrr_0.3.4     readr_2.1.2     tidyr_1.2.0     tibble_3.1.8   
## [13] ggplot2_3.3.6   tidyverse_1.3.2 plyr_1.8.7     
## 
## loaded via a namespace (and not attached):
##  [1] TH.data_1.1-1       googledrive_2.0.0   minqa_1.2.4        
##  [4] colorspace_2.0-3    ggsignif_0.6.3      ellipsis_0.3.2     
##  [7] estimability_1.4    base64enc_0.1-3     fs_1.5.2           
## [10] rstudioapi_0.13     farver_2.1.1        fansi_1.0.3        
## [13] mvtnorm_1.1-3       lubridate_1.8.0     xml2_1.3.3         
## [16] codetools_0.2-18    splines_4.2.1       cachem_1.0.6       
## [19] knitr_1.39          jsonlite_1.8.0      nloptr_2.0.3       
## [22] pbkrtest_0.5.1      broom_1.0.0         dbplyr_2.2.1       
## [25] compiler_4.2.1      httr_1.4.3          emmeans_1.7.5      
## [28] backports_1.4.1     assertthat_0.2.1    fastmap_1.1.0      
## [31] gargle_1.2.0        cli_3.3.0           htmltools_0.5.3    
## [34] tools_4.2.1         coda_0.19-4         gtable_0.3.0       
## [37] glue_1.6.2          Rcpp_1.0.9          carData_3.0-5      
## [40] cellranger_1.1.0    jquerylib_0.1.4     vctrs_0.4.1        
## [43] nlme_3.1-157        xfun_0.31           rvest_1.0.2        
## [46] lifecycle_1.0.1     rstatix_0.7.0       googlesheets4_1.0.0
## [49] MASS_7.3-57         zoo_1.8-10          scales_1.2.0       
## [52] hms_1.1.1           parallel_4.2.1      sandwich_3.0-2     
## [55] RColorBrewer_1.1-3  yaml_2.3.5          gridExtra_2.3      
## [58] gdtools_0.2.4       sass_0.4.2          stringi_1.7.8      
## [61] highr_0.9           zip_2.2.0           boot_1.3-28        
## [64] systemfonts_1.0.4   rlang_1.0.4         pkgconfig_2.0.3    
## [67] evaluate_0.15       lattice_0.20-45     labeling_0.4.2     
## [70] cowplot_1.1.1       tidyselect_1.1.2    magrittr_2.0.3     
## [73] R6_2.5.1            generics_0.1.3      multcomp_1.4-19    
## [76] DBI_1.1.3           pillar_1.8.0        haven_2.5.0        
## [79] withr_2.5.0         mgcv_1.8-40         survival_3.3-1     
## [82] abind_1.4-5         modelr_0.1.8        crayon_1.5.1       
## [85] car_3.1-0           uuid_1.1-0          utf8_1.2.2         
## [88] officer_0.4.3       tzdb_0.3.0          rmarkdown_2.14     
## [91] grid_4.2.1          readxl_1.4.0        data.table_1.14.2  
## [94] reprex_2.0.1        digest_0.6.29       xtable_1.8-4       
## [97] numDeriv_2016.8-1.1 munsell_0.5.0       bslib_0.4.0
rm(list=ls())